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Abstract 

f 

This  program  has  concentrated  on  three  major  areas:  the  application  of  high 
angular  resolution  image  reconstruction  techniques  to  the  reconstruction  of  solar  surface 
features;  speckle  imaging  of  a  wide  range  of  astronomical  sources;  and  the  implementation 
of  adaptive  optics  for  faint  object  imaging.  In  this  project,  we  have  made  substantial 
modifications  to  a  technique  for  reconstruction  of  high  resolution  images  from  single  short 
exposure  solar  frames,  blind  iterative  deconvolution  (IDC).  We  have  been  testing  these 
improvements  using  numerical  simulation  data.  Analysis  of  speckle  data  of  the  supernova 
SN1987A  has  detected  a  new  bright  source,  0.9  arcseconds  south  of  the  SN,  as  well  as 
substantial  structure  in  the  region  surrounding  the  SN.  We  also  have  new  results  on 
several  Young  Stellar  Objects  and  supergiants.  In  the  area  of  adaptive  optics,  an  AOA 
wavefront  sensor  has  been  set  up  and  tested  using  an  image  inten  sifter  which  increases  its 
sensitivity  by  three  orders  of  magnitude.  Faint  object  image  active  tilt  correction  has  also 
been  tested  with  some  important  enhancements,  including  new  highly  linear  and  sensitive 
CCD  quad  cells  developed  by  Cal  Tech  and  Tektronix,  and  an  off-the-shelf  high  speed  2-D 
tilting  mirror  with  greatly  improved  specifications.  The  effects  of  only  partially  correcting 
atmospheric  turbulence  have  been  theoretically  analyzed  and  numerically  simulated.  [  ' 

< 

I.  SOLAR  IMAGING 


Extensive  experiments  performed  in  the  area  of  solar  speckle  imaging  have  proven 
to  be  limited  in  the  accuracy  of  the  reconstructions  due  to  the  long  sequences  of  data 
required  for  diffraction  limited  recovery.  Solar  surface  features  appear  to  evolve  on  time 
scales  substantially  shorter  than  the  time  required  to  take  a  long  enough  series  of  frames 
with  a  CCD  camera.  This  problem  was  demonstrated  through  experiments  which  tested 
reconstruction  algorithms  using  computer  generated  simulation  data,  and  by  processing 
shorter  sets  of  real  solar  data,  "lhe  non-time-varying  simulations  '^.covered  the  highest 
frequency  information  with  excellent  accuracy.  In  the  reJ  short  data  sets  (20  frames 
@  2  frames/sec)  recovered  artifact  free  but  unconvergcd  images,  while  longer  data  sets 
(100  frames)  produced  fringe  like  artifacts  in  the  reconstructions,  apparently  due  to 
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evolution  in  the  solar  features.  This  rapid  evolution  of  high  resolution  features  has  recently 
been  documented  in  Skylab  solar  data  sequences,  and  in  sequences  of  solar  frames 
recorded  under  optimal  seeing  conditions  at  the  La  Palma  solar  telescope. 

Obviously,  one  possible  approach  to  solving  these  problems  is  a  more  rapid 
recording  of  data.  However,  obtaining  long  enough  sequences  under  good  conditions  is 
difficult  and  could  not  be  used  in  any  sort  of  routine  observing  mode.  In  this  program,  we 
have  directed  out  efforts  toward  the  development  of  algorithms  which  could  reconstruct 
images  from  much  shorter  sequences  of  data.  In  solar  imaging,  the  signal-to-noise  of  each 
frame  is  very  high,  due  to  the  high  light  levels,  and  this  is  the  major  requirement  for  the 
successful  application  of  Blind  Iterative  Deconvolution  techniques. 

Blind  Iterative  Deconvolution 

Blind  Iterative  deconvolution  (ID)  (Ayers  and  Dainty,  1988)  combines  iterative 
techniques  such  as  those  developed  for  phase  retrieval  (Fienup,  1978)  with  blind 
deconvolution  (Lane  and  Bates,  1987).  One  starts  with  an  image  which  is  degraded  by 
some  blurring  function.  A  necessary  condition  for  the  algorithm  to  work  is  that  the 
blurring  function  be  constant  over  the  entire  image  field  to  be  restored  (stationarity).  It  is 
also  assumed  that  the  degradation  is  a  linear  operation.  The  general  approach  is  then  to  find 
a  pair  of  functions  whose  convolution  gives  the  input  image  within  a  set  of  physical 
constraints,  such  as  positivity  of  the  two  convolved  functions  or  the  signal-to-noisc  ratio  in 
the  Fourier  transform  (FT).  While  it  has  not  been  proven  that  the  derived  functions  are 
unique,  complicated  images  appear  to  converge  on  only  one  sensible  solution. 
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A  flow  diagram  for  the  technique  is  given  in  figure  1 .  One  starts  with  a  degraded 
image  and  an  initial  estimate  of  the  point  spread  function  (PSF).  The  initial  PSF  can  be 
completely  random,  however  the  number  of  cycles  of  iteration  required  to  converge  on  an 
acceptable  answer  is  highly  dependant  on  how  close  the  first  estimate  of  the  PSF  is  to  the 
actual  PSF.  Both  inputs  are  Fourier  transformed  and  a  deconvolution  is  performed  by 
dividing  the  FT  of  the  image  by  the  FT  of  the  PSF.  Constraints  are  applied  to  the  ratio 
based  on  estimated  signal-to-noise  ratios.  The  result  is  transformed  back  to  image  space 
and  the  positivity  and  support  constraints  are  applied  (the  image  support  is  defined  as  the 
region  in  which  the  image  is  non-zero).  The  negatives  in  the  image  are  summed  and  the 
pixels  with  negative  values  are  set  to  zero.  The  negatives  are  then  uniformly  subtracted 
within  the  support  region  in  order  to  preserve  the  total  power  in  the  image.  The  sum  of  the 
negatives  is  also  used  as  a  diagnostic  of  convergence.  The  FT  of  the  original  degraded 
image  is  then  deconvolved  by  the  FT  of  the  image  obtained  from  the  previous  iteration. 
Fourier  plane  constraints  are  applied  and  the  result  is  transformed  back  to  image  space. 
Again,  positivity  and  the  support  constraint  are  enforced.  The  result  is  a  new  estimate  of 
the  PSF.  The  iteration  continues  until  some  concluding  criterion  is  reached.  Either  a  limit 
to  the  number  of  cycles,  or  the  level  of  residual  negatives  in  the  image,  or  some  other  image 
quality  criterion  may  determine  the  stopping  point  After  examination  of  the  output  image 
and  PSF,  the  results  may  be  fed  back  into  the  loop  for  continued  iterations. 

There  are  a  number  of  parameters  which  must  be  chosen  in  order  to  ensure 
convergence  and  an  optimum  result.  Probably  the  most  important  are  estimates  of  the 
signal-to-noisc  ratio  in  the  data  and  the  range  of  spatial  frequencies  over  which  the  image 
recovery  operation  is  performed.  A  low  pass  filter  is  applied  in  frequency  space  after  each 
cycle  and  then  its  inverse  filter  is  applied  before  each  deconvolution.  The  extent  of  this 
filter  is  adjusted  for  the  signal-to-noise  ratio  in  the  image  and  the  PSF. 
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Reconstructions  with  ID 

The  work  on  ID  already  undertaken  at  CfA  has  produced  an  algorithm  that  has  been 
tested  with  computer  simulations  and  also  applied  to  some  real  data.  Results  of  the 
simulations  are  shown  in  Figs.  2  and  3.  Fig.  2a  shows  the  input  diffraction  limited  image 
of  8  point  sources,  the  bottom  right-hand  "point"  being  two  unresolved  points.  This  image 
was  convolved  with  the  PSF  in  Fig.  2b  to  produce  the  input  image  to  be  deconvolved  with 
ED  shown  in  Fig.  2c.  The  starting  guess  for  the  PSF  (Fig.  2d)  used  to  perform  the  first 
cycle  of  ID  was  a  gaussian  with  random  noise  and  a  half  power  width  with  approximately 
the  same  radius  as  the  input  image.  In  most  real  situations,  there  is  usually  some 
reasonable  estimate  of  the  PSF  which,  when  used  as  a  first  guess,  should  improve  the  rate 
of  convergence.  Here,  a  random  PSF  was  used  to  demonstrate  the  dramatic  evolution  of 
the  reconstructed  PSF  towards  the  actual  diffraction  limited  PSF,  despite  the  poor  starting 
estimate.  The  image  and  PSF  obtained  from  ID  after  20  cycles  are  shown  in  Figs.  3a  and 
3b,  respectively.  The  major  features  present  in  the  diffraction  limited  image  and  PSF  have 
been  reconstructed,  however  further  iterations  lead  to  a  considerable  improvement  in  the 
reconstructions.  Figs.  3c  and  3d  show  the  image  and  PSF  reconstructions  after  180  cycles 
of  ID.  Comparison  of  these  results  with  Figs.  2a  and  2b  show  a  dramatic  recovery  of  both 
the  morphology  and  intensities  present  in  the  diffraction  limited  image  and  PSF. 

These  initial  computer  simulations  have  produced  highly  encouraging  results. 
However,  detailed  characterization  of  the  effects  of  noise  and  the  robustness  of  the 
convergence  under  variations  in  reconstruction  parameters  have  not  been  rigorously 
investigated.  For  ID  or  any  deconvolution  method  to  be  a  useful  tool,  it  must  be 
rigorously  characterized  with  both  simulations  and  real  data. 


Some  initial  attempts  have  been  made  to  use  the  technique  on  real  x-ray  and  optical 
CCD  data.  Figs.  4a-c  show  the  results  of  processing  a  128x128  pixel  Einstein  image  of 
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SNR  0540-69.3.  Fig.  4a  is  a  contour  plot  of  the  input  image  (the  input  image  had  one 
smoothing  pass).  Fig.  4b  shows  the  result  from  ID  after  only  10  cycles  of  iteration.  There 
has  obviously  been  a  great  increase  in  the  intensity  and  sharpening  of  the  central  source. 
Structure  around  the  source  is  evident  at  scales  of  about  10  arcseconds.  Fig.  4c  shows  the 
recovered  PSF.  Additional  cycles  of  ID  made  almost  no  change  in  the  reconstruction, 
indicating  a  stable  solution  for  the  chosen  parameters.  It’s  possible  that  other  parameters, 
such  as  allowing  the  signal-to-noise  ratio  a  greater  latitude  might  produce  some  additional 
sharpening.  However  this  result  is  already  very  promising.  It  will  be  interesting  to 
compare  these  results  to  optical  images  processed  in  a  similar  manner,  where  extended 
structure  is  known  to  be  on  the  same  scale  as  the  recovered  features. 

Figs.  5a-c  show  results  of  ID  processing  on  CCD  images  of  supernova  SN1987A 
recorded  at  the  CTIO  4-mcter  telescope  in  March,  1989  (Heathcote,  1989).  5a  shows  the 
original  image  of  the  supernova,  its  two  companions,  and  the  suspected  optical  echo.  The 
wide  companion  is  separated  by  2.9  arcseconds  from  the  SN.  5b  is  the  ID  reconstruction, 
again  after  only  10  cycles  of  processing  and  5c  shows  the  derived  PSF.  Clearly  the  image 
has  been  sharpened  and  both  the  stars  and  the  echo  are  now  well  resolved.  These  results 
demonstrate  the  power  of  the  technique,  and  show  that  one  obtains  rapid  convergence  if  a 
good  initial  estimate  of  PSF  is  available. 

We  are  currendy  testing  approaches  to  adding  what  may  be  a  still  more  powerful 
convergence  constraint.  If  a  series  of  short  exposure  images  is  recorded  where  the  object 
remains  constant  but  the  degrading  atmosphere  changes,  then  each  frame  should  converge 
to  the  same  image,  but  with  very  different  PSFs.  Forcing  commonalty  between  separate 
results  from  two  or  more  independent  iterations  should  prove  to  have  an  important  effect  on 
convergence  properties. 
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Fig.  3c  The  reconstruction  of  the  image  after  180  cycle*  of 
iterative  deconvolution. 


Fig.  3d  The  reconstruction  of  the  PSF  after  180  cycles  ol 
iterative  deconvolution. 
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In  the  solar  imaging  problem,  we  had  only  limited  success  using  standard  speckle 
image  reconstruction  techniques  for  reconstructing  high  resolution  solar  images,  due  to  the 
large  number  of  frames  required  to  give  good  signal-to-noise  ratio  images.  With  the  IDC 
technique,  only  a  few  frames  are  required,  so  the  stability  problem  is  not  an  issue.  Other 
questions  that  remain  are  the  problem  of  isoplanicity  (over  how  large  an  angle  is  the 
atmospheric  PSF  stationary?),  and  do  the  boundaries  of  the  image  caused  by  the  Finite  size 
of  the  detector  cause  degradation?  Some  of  the  techniques  we  have  used  for  speckle  will 
play  an  important  role,  particularly  "detrending  ',  where  edge  effects  are  reduced  by 
subtraction  of  a  low-pass  filtered  version  of  the  image  before  processing.  Also,  what  are 
the  effects  of  using  an  imperfect  detector?  The  isoplanatic  problem  and  the  detector  effects 
can  only  be  tested  with  real  data,  though  we  are  able  to  add  simple  forms  of  noise  to  the 
simulations,  while  the  boundary  problem  may  be  explored  with  simulations. 

Tests  of  the  IDC  approach  arc  continuing  with  CCD  data  recorded  at  the 
Sacramento  Peak  Solar  Observatory  and  data  recorded  by  a  group  from  Lockheed  who 
observed  using  the  Swedish  solar  telescope  at  La  Palma  .  The  atmospheric  conditions  are 
sufficiently  good  so  that  some  individual  frames  are  nearly  diffraction  limited.  By 
obtaining  sequences  of  frames  where  the  seeing  is  very  good  and  then  degrades,  one  can 
compare  the  results  to  the  very  good  frames,  giving  a  form  of  "ground  troth"  for  testing 
the  IDC  reconstruction  process.  Processing  several  separate  degraded  frames  and 
comparing  the  results  to  each  other  and  the  best  frames  is  also  important 

We  believe  that  ID  could  produce  very  good  solar  physics.  It  could  allow  long 
sequences  of  frames  to  be  corrected,  and  mosaics  of  much  larger  fields  than  the  isoplanatic 
patch  could  be  synthesized  from  overlapping  reconstructions  of  small  areas.  It  could  also 
be  applied  to  data  recorded  at  telescopes  with  large  apertures,  where  the  atmosphere  will 
always  have  a  substantial  degrading  effect 
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ll.  ADAPTIVE  OPTICS 

Star-Tracker 

A  new  adaptive  optics  laboratory  has  been  equipped  and  is  now  in  operation. 
Initially  it  was  used  to  evaluate  the  performance  of  an  experimental  "star-tracker”,  which 
consists  of  a  position  sensor,  a  processor  and  a  simple  actuator.  These  components  form  a 
closed  loop  control  system.  The  position  of  an  attenuated  laser  spot  is  monitored  by  a  quad¬ 
cell  and  the  resulting  four  voltages  are  read  by  a  Macintosh  II  computer.  A  Kalman  filter 
has  been  implemented  on  the  Macintosh  n  which  predicts  the  position  of  the  laser  spot 
from  its  previous  displacements  across  the  quad-cell.  From  the  output  of  this  Kalman  filter, 
a  control  voltage  (ranging  from  0  to  +10  volts)  is  derived  which  is  then  used  to  control  the 
piezoelectric  actuators  and  return  the  laser  spot  to  the  center  of  the  quad-cell.  The  position 
of  the  laser  spot  is  controlled  by  a  gimbal  mounted  mirror  which  can  be  tilted  about  two 
perpendicular  axes  by  two  piezostatic  actuators.  However,  as  the  actuators  produce  their 
maximum  displacement  with  a  driving  voltage  of  1000  volts  an  amplifier  has  been 
constructed  to  multiply  the  control  voltage  by  a  factor  of  one  hundred  before  it  is  applied  to 
the  actuators.  The  Mac  II  is  also  used  to  produce  atmospheric  simulation  error  signals 
which  can  be  fed  to  a  second  active  minor  in  the  system.  These  errors  are  independent  and 
random,  with  statistics  matching  an  atmospheric  model.  The  correction  system  then 
performs  its  operation  and  the  residuals  are  measured.  This  allows  a  full  test  and 
optimization  of  the  tracking  system  in  the  lab,  over  a  realistic  set  of  conditions,  at  a  range  of 
light  levels,  time  constants,  and  atmospheric  variances. 

Two  new  hardware  components  should  aid  in  the  implementation  of  faint  object 
high  speed  active  guiding.  Special  CCD  quad  cells  have  been  developed  by  a  collaboration 
between  groups  at  Tektronix,  Cal  Tech,  and  Johns  Hopkins  University.  These  cells  have 
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very  high  quantum  efficiency,  very  low  readout  noise  (a  few  electrons)  and  can  be  read  out 
at  kilohertz  rates.  They  appear  to  be  ideal  for  detection  of  image  centroids  at  the  low  light 
levels,  required  for  astronomical  rapid  guiding  and  interferometry.  Cal  Tech  is  making 
some  of  these  chips  available  to  us,  and  Johns  Hopkins  is  providing  designs  for  mounting 
and  cooling  the  chips.  We  have  also  located  a  high  speed  tilt  mirror  which  should  have 
greatly  improved  performance  compared  to  the  Burleigh  minors  which  we  have  been  trying 
;o  use  for  this  purpose.  The  mirror  was  designed  by  a  German  company,  Physike  Corp., 
for  use  by  the  European  Southern  Observatory.  There  are  a  few  specifications  that  we  are 
still  checking  on  (vibration  and  piston  motion)  but  from  what  we  have  already  learned,  this 
mirror  appears  to  have  more  than  adequate  characteristics  for  high  speed  image  tracking  and 
correction. 

AOA  Wavefront  Sensor 

A  wavefront  sensor,  which  uses  the  Hartmann  test  to  determine  atmospheric 
wavefront  distortions,  has  been  constructed  and  demonstrated  by  Adaptive  Optics 
Associates  (AOA)  under  the  control  of  their  MassComp  computer.  AOA  then  undertook  to 
port  the  controlling  software  from  MassComp  Unix  to  VMS  on  our  Dec  MicroVax. 
Unfortunately  their  experience  with  VMS  was  very  limited  and  they  were  over  optimistic  in 
the  assessment  of  their  ability  to  translate  the  operational  environment  from  Unix  to  VMS. 
Therefore,  it  was  decided  to  install  a  Unix  operating  system  on  the  MicroVax.  Once  this 
was  completed  the  graphics,  data  transfer,  and  control  software  were  installed  with  little 
difficulty  and  A.O.A.  The  wavefront  sensor  has  now  been  delivered  and  initial  tests  have 
been  performed  with  it  The  next  step  will  be  to  produce  a  set  of  controlled  aberrated 
wavefronts  to  test  our  approach  of  computational  estimation  and  correction.  We  will  then 
extend  these  test  by  interfacing  the  sensor  to  our  atmospheric  simu’ator  and  photon  camera. 
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This  will  allow  a  full  test  of  our  partial  atmospheric  correction  and  adaptive  optics-  speckle 
imaging  applications. 

We  have  completed  the  implementation  of  improved  software  on  the  Microvax 
which  controls  the  AOA  wavefront  sensor.  The  use  of  correlation  techniques  to  locate  and 
apply  the  array  of  subapertures  has  proven  to  be  very  effective.  With  the  software,  as 
delivered  by  AOA,  determination  of  the  positions  of  the  subapertures  was  extremely  erratic. 
Since  this  is  a  critical  element  in  the  operation  of  the  sensor,  testing  could  not  be  continued 
until  this  problem  was  solved.  There  were  numerous  other  problems  with  the  control  and 
display  software,  some  of  them  due  to  the  conversion  of  the  software  from  a  Masscomp 
computer  to  the  Vax,  and  some  fundamental  to  the  software  design.  The  software  is  finally 
in  sufficiendy  good  shape  to  allow  more  extensive  laboratory  tests  of  faint  object  wavefront 
sensing  and  correction.  Problems  with  the  three  stage  image  intensifier ,  needed  to  allow 
the  sensor  to  work  at  photon  counting  levels  for  astronomical  sources,  were  finally 
corrected  during  this  reporting  period.  Extensive  tests  of  the  low  light  mode  of  operation 
have  not  yet  been  performed. 

Partial  Atmospheric  Correction  Analysis 

Wavefront  sensors  such  as  the  AOA  system  we  have  purchased  arc  expensive  and 
technologically  complicated,  requiring  massive  parallel  computation  and  extensive 
electronic  control  systems  for  their  operation.  A  large  telescope  under  typical  seeing 
conditions  will  have  several  hundred  to  thousands  of  correlation  areas  across  its  pupil, 
requiring  an  equal  number  of  correction  channels  in  an  adaptive  optics  system  if  full 
correction  is  to  be  obtained.  Our  system  has  only  37  correction  elements,  so  it  could  only 
fully  correct  a  very  small  telescope.  However,  we  have  performed  numerical  experiments 
and  analyses  to  determine  the  effectiveness  of  using  a  corrector  with  many  fewer  elements 


-13- 


than  the  number  of  correlation  areas  in  the  pupil.  The  results  show  (Nisenson  and  Barakat, 
1987)  that  even  a  very  limited  number  of  elements  can  produce  a  surprisingly  large 
improvement  in  the  atmospheric  "seeing".  The  wavefront  sensor,  in  effect,  corrects  only 
the  lower  spatial  frequencies,  leaving  a  residual  high  frequency  background,  but 
concentrating  the  energy  in  the  system  point  spread  function.  Correcting  only  low  spatial 
frequencies  allows  light  collection  over  larger  regions  in  the  pupil  and  these  regions  have  a 
commensurately  longer  time  constant.  This  allow  the  system  to  operate  to  substantially 
fainter  magnitude  limits.  The  detailed  results  of  this  analysis  are  included  as  an  appendix  in 
this  report. 


III.  SPECKLE  IMAGING  ASTRONOMY 

The  New  PAPA  2-d  Photon  Counting  Detector 

A  new  version  of  the  PAPA  2-dimensional  photon  counting  camera  has  been  developed,  in 
collaboration  with  four  groups  at  other  institutions  who  are  involved  with  speckle  imaging  and 
long-baseline  interferometry.  Five  cameras  have  been  built  with  greatly  improved  characteristics 
compared  to  the  original  PAPA,  which  was  built  with  funding  from  AFOSR.  The  primary 
advantages  of  the  new  PAPA  are  increased  quantum  efficiency  (>10  %),  higher  maximum  count 
rate  (>1  Mhz),  a  much  lower  error  rate,  and  much  better  thermal  and  mechanical  stability.  The  new 
camera  also  has  1 1  bit  by  1 1  bit  resolution  (2048x2048  pixels)  compared  to  the  256x256  resolution 
of  the  old  camera.  The  increased  sensitivity  and  speed  of  the  new  camera  allows  us  to  obtain 
enhanced  resolution  imaging  of  the  SN  and  its  circumstellar  environment,  despite  their  overall 
decline  in  brightness. 
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Observations  of  SN1987A 

High  angular  resolution  observations  of  SN1987A  (SN)  have  provided  valuable 
data  that  contribute  to  the  understanding  of  this  extraordinary  astrophysical  event.  Using 
speckle  interferometric  techniques,  it  is  possible  to  overcome  the  resolution  limits  imposed 
by  atmospheric  turbulence  and  observe,  in  some  detail,  the  morphology  of  the  SN  and  its 
environment.  A  wide  variety  of  other  techniques  such  as  spectroscopy,  photometry, 
polarimetry,  X-ray  and  gamma-ray  observations  are  providing  an  extensive  data  base  that, 
combined  with  the  high  angular  resolution  measurements  will  give  a  detailed  understanding 
of  the  mechanisms  and  dynamics  of  supemovae,  testing  the  validity  of  theoretical  models. 

Speckle  observations  at  CTIO  in  March- April  revealed  the  presence  of  a  second 
bright  source  separated  by  only  60  mas  from  the  SN,  and  within  a  factor  of  10  of  its 
brightness.  Furthermore,  analysis  of  data  acquired  in  March-April  along  with  data  from 
two  other  runs  in  May-June,  and  November,  1987,  provided  measurements  of  the  apparent 
diameter  of  the  SN  at  several  wavelengths.  Most  recently,  data  recorded  in  March  and 
April  of  1988  show  evidence  for  strong  asymmetries  in  the  expanding  SN  shell,  with  the 
major  axis  of  the  asymmetry  in  the  same  direction  as  the  measured  axis  of  polarization,  as 
well  as  the  direction  of  earlier  detected  bright  2nd  source. 

Diameter  measurements  were  made  by  fitting  the  integrated  power  spectra  obtained 
from  the  speckle  process  to  the  power  spectrum  of  a  uniform  disk  (an  Airy  disk).  Accurate 
fitting  allows  diameter  estimates  at  scales  well  below  the  "diffraction  limit"  of  the  telescope, 
to  accuracies  of  a  few  mi  Ui  arc  seconds. 

While  measurements  of  equivalent  precision  have  not  previously  been  demonstrated 
using  speckle  techniques,  the  validity  of  our  results  has  been  supported  by  laboratory 
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simulations  and  measurements  of  stars  whose  diameters  have  been  determined  by  other 
techniques.  Results  show  that  diameters  measured  in  Ha  are  in  good  agreement  with 
predictions  based  on  measurements  of  the  velocity  from  the  Ha  trough.  Measurements  in 
the  earliest  months  after  the  explosion  showed  a  surprisingly  large  diameter  in  the 
continuum,  but  results  in  the  later  months  gave  results  results  closer  to  predicted  scales. 
We  have  set  up  an  extended  observing  program  in  collaboration  with  scientists  at  CTIO,  so 
that  speckle  data  is  being  acquired  on  once-a-month  basis.  Latest  results  indicate  a 
substantial  asymmetry  in  the  supernova  shell,  with  an  axis  aligned  in  the  direction  of  the 
measuredpolarization.  As  the  supernova  expands,  its  details  should  be  observable  with  an 
ever  increasing  degree  of  detail.  Our  technique  appears  to  be  providing  a  unique  window 
on  this  one-of-a-kind  object.  The  success  of  our  measurements  to  better  than  1  mas 
precision  may  also  provide  a  tool  for  many  other  scientific  problems  which  was  previously 
thought  only  possible  with  long  baseline  interferometers,  but  at  much  fainter  magnitude 
limits  than  would  be  achievable  with  an  interferometer. 

Reduction  of  speckle  data  recorded  on  November  28,  1988,  December  20,  1988, 
and  January  20, 1989  show  the  existence  of  a  bright  source  near  SN1987A  (LAU  Circulars 
#  4749  and  4752).  Observations  performed  with  the  CTIO  4-meter  telescope,  the  PAPA 
detector,  and  speckle  image  reconstruction  techniques  show  a  source  located  at  P.A.  200 
(+/-  5)  degrees  for  all  observations.  Preliminary  estimates  of  the  magnitude  differences 
from  the  SN  (551  nm)  and  angular  separations  were:  3.3  (+/-  0.2)  magnitudes  and  0.85 
(+/-  0.05)  arcsec  (November  28);  2.6  (+/-  0.2)  and  0.85  (+/-  0.05)  (December  20);  and  2.1 
(+/-  0.2)  magnitudes  and  0.95  (+/-  0.05)  arcsec  (January  20).  The  separation  and  position 
angle  of  this  source  is  consistent  with  the  position  of  the  compact  blob  reported  by  Allen  et 
al  (LAU  circular  #  4747).  It  is  also  consistent  with  a  position  predicted  for  uniform  motion 
of  the  bright  source  which  was  detected  using  speckle  techniques  on  April  1,  1987  (1AU 
circular#  4382  )  and  by  Matcher  et  al  on  April  15,  1987  (IAU  Circular#  4413).  Wc  note 
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that  the  brightness  of  this  source  is  consistent  with  the  rate  of  decay  of  the  supernova,  so  if 
this  is  the  same  object  detected  almost  two  years  ago,  it's  energy  source  may  well  be 
Cobalt-56  decay. 

We  were  also  able  to  detect  star  2  separated  by  2.8  arc  seconds  from  the  SN  with  a 
magnitude  difference  of  about  3  (at  551  nm)  and  position  angle  of  315  (+/-  5)  degrees  and 
star  3  is  just  barely  detectable  at  a  separation  of  1.7  arcseconds  and  position  angle  1 10 
degrees.  We  did  not  detect  any  feature  to  the  north  of  the  SN,  as  suggested  by  Allen’s 
long-slit  spectral  measurements.  Initial  reductions  of  CCD  data  recorded  at  the  prime  focus 
of  the  CTIO  4-meter  show  possible  confirmation  of  these  detections. 

It  should  continue  to  be  extremely  interesting  to  observe  the  SN  as  it  evolves,  in 
particular  when  it  allows  a  first  close  look  at  the  formation  of  an  SNR.  We  will  look  for 
scattering  sources  and  possible  evidence  for  early  signs  of  synchrotron  radiation.  If  the 
rapidly  spinning  pulsar  is  actually  in  the  center  of  the  SN,  as  early  observations  suggest, 
than  speckle-polarization  experiments  may  be  able  to  detect  its  effect  on  the  surrounding 
medium.  Finally,  SN1987A  has  provided  many  unexpected  results  during  the  past  3  years 
with  behavior  which  has  not  been  predicted  by  existing  SN  models.  We  believe  that  this  is 
a  compelling  reason  for  continued  careful  monitoring  with  all  available  techniques. 

Young  Stellar  Objects 

Results  from  our  study  of  several  YSO’s  show  a  number  of  stars  with  extended 
disk -like  structures.  Observations  of  DG  Tau  at  several  wavelengths  including  6308  A 
[01],  the  adjacent  continuum  and  Ha  reveal  extended  asymmetric  structure  at  angular  scales 
larger  than  0.3  arcsecond.  This  asymmetry  is  located  along  NW-SE  direction  and  is 
aligned  with  the  direction  of  the  observed  jet.  An  extended  halo  around  Walker  90  was 
also  discovered.  The  halo  is  elongated  in  a  SE-NW  direction  and  is  almost  perpendicular  to 
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the  measured  polarization  angle.  Results  obtained  from  fitting  a  power  law  model  to  the 
intensity  profile  of  the  extended  halo  indicates  the  presence  of  a  flared  disk  The 
obscuration  of  the  central  star  was  also  obtained  from  the  fit.  S  CrA  is  a  close  (one 
arcsecond)  separation  binary  in  which  both  components  are  pre-main  sequence  stars.  Our 
observations  reveal  extended  disk-like  structures  surrounding  both  components  of  the 
binary  system.  There  also  appears  to  be  either  a  knot  of  material  between  the  two  stars  or  a 
faint  third  component  in  the  system. 

Late-type  Supergiants 

Several  new  measurements  of  the  position  angle  of  the  close  companion  to  the 
supergiant  a  Ori  have  been  made.  We  have  also  discovered  a  large  bright  structure  (or  few 
unresolved  structures)  on  the  surface  of  a  Ori  which  are  probably  related  to  the 
phenomenon  of  stellar  super-granulation. 

We  have  detected  an  elongation  of  o  Ceti's  (Mira)  upper  atmosphere  (from 
observations  obtained  at  two  phases  of  the  pulsation  cycle  of  the  star).  This  structure  is 
probably  associated  with  a  disk  around  the  star  or  a  molecular  cloud  in  the  vicinity  of  the 
star  formed  as  a  result  of  stellar  pulsations. 
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Adaptive  optics  has  9een  only  limited  application  in  astronomical  facilities,  despite  its  significant  potential  for 
improving  seeing  conditions  and  increasing  observing  efficiency  and  prodnrti'ip1  Expense  ami  technological 
difficulty  appear  to  be  the  reasons  that  this  is  the  case.  Correction  of  large  apertures  requires  hundreds  of  active 
elements  in  both  the  wave  front  sensor  and  the  adaptive  mirror.  We  have  performed  some  one  dimensional 
numerical  simulations  to  test  atmospheric  wave  front  correction  when  the  active  element  is  not  matched  to  the 
correlation  scale  in  the  pupil  The  results  demonstrate  that  substantial  seeing  improvement  can  be  obtained  with 
an  adaptive  optical  system  having  a  limited  number  of  active  elements. 


1.  INTRODUCTION 

The  technology  for  the  real-time  correction  of  atmospheric 
aberrations  using  adaptive  optics  has  been  under  develop¬ 
ment  for  more  than  a  decade. However,  except  for  a  few 
limited  cases,  these  techniques  have  not  been  applied  to 
uptical  ostrunomy,  probably  because  of  the  difficulty  and 
expense  of  the  necessary  technologies.  Still,  a  technique 
that  improves  the  atmospheric  seeing  would  be  of  great 
benefit  for  almost  all  astronomical  scientific  programs,  in¬ 
creasing  magnitude  limits,  signal-to-noise  ratio,  and  the  ob¬ 
servational  efficiency. 

Virtually  all  adaptive  optical  systems  considered  for  see¬ 
ing  correction  have  been  analyzed  for  the  case  when  the  scale 
of  the  active  element  (in  the  wave  front  sensor  and  the 
adaptive  mirror)  matches  the  scale  of  the  atmospheric  corre¬ 
lation  length  (r0).  For  large  telescopes,  this  means  that, 
even  for  very  good  seeing  (i.e-,  <1  arcstc),  hundreds  of  cor¬ 
rection  elements  are  needed.5  For  example,  1-arcsec  seeing 
with  a  2-m  telescope  would  require  more  than  400  active 
elements  for  diffraction-limited  performance.  When  one 
combines  this  large  number  of  active  elements  with  the  re¬ 
quirement  for  high-speed  operation  (1-msec  cycle  time),  the 
adaptive  optical  system  becomes  prohibitively  expensive. 
In  addition,  even  an  optimized  system  will  have  a  magnitude 
limit  that  restricts  its  utility  to  bright  sources  and  a  limited 
class  of  scientific  problems. 

An  alternative  approach  is  to  consider  the  effects  of  an 
adaptive  optical  system  with  active  elements  having  a  scale 
that  is  much  larger  than  r0,  with  the  goal  of  obtaining  partial 
improvement  in  the  seeing.  For  many  astronomical  prob¬ 
lems  even  a  limited  improvement  in  seeing  would  yield  sub- 
a.jntial  gains,  reducing  data  integration  times  and  enhanc¬ 
ing  image  quality.  Since  each  active  element  is  increased  in 
area,  many  more  photons  are  detected  per  element,  improv¬ 
ing  the  magnitude  limit  of  the  wave-front  measurement  pro¬ 
portionately.  The  longer  scales  in  the  atmosphere  should 
also  have  longer  time  constants,  so  the  system's  cycle  time 


would  increase,  again  improving  the  magnitude  limit  and 
reducing  the  technological  requirements.  Finally,  correc¬ 
tion  of  only  the  longer-seale  atmospheric  structure  could 
also  result  in  a  wider  isoplanatic  field,  increasing  the  proba¬ 
bility  of  finding  a  sufficiently  bright  offset  star  for  wave- 
front  measurement,  permitting  correction  and  integration 
on  faint  or  extended  sources. 

In  this  paper  we  describe  some  numerical  experiments 
that  we  have  performed  to  test  the  concept  of  partial  correc¬ 
tion.  We  have  generated  some  simulated  one-dimensional 
atmospheric  wave  fronts,  applied  correction  factors  to  the 
wave  fronts  for  various  spatial  Beales,  and  then  calculated 
transfer  functions,  point-spread  functions  (psfs),  and  im¬ 
ages  for  intercomparison.  The  results  support  the  thesis 
that  substantial  improvement  in  imaging  capabilities  could 
be  obtained,  even  with  limited  adaptive  atmospheric  correc¬ 
tion. 


2.  WAVE-FRONT  GENERATION  AND 
CORRECTION 

The  psf6  for  a  slit  aperture,  without  scintillation  or  deter¬ 
ministic  aberrations,  is 


Hr)  = 


1  f*' 

—  exp[i)tZ(p)]exp(u'p)dp 

2  l-i 


(2.1) 


Zip)  is  the  spatially  random  wave-front  aberration  function 
over  the  exit  pupil  induced  by  the  turbulent  atmosphere. 
We  model  Zip)  as  a  zero-mean,  spatially  atationarv  random 

process; 

(Zip))  =  0,  (2.2) 

<Z(p,)Z(pj)>  =  o2r(|p2  -  p,|),  (2.3) 

where  a2  is  the  variance  of  the  random  wave  front  measured 
in  X  units  and  r(Ap)  is  the  wave-front  correlation  function 
(the  generation  of  sample  realizations  of  Zip )  is  outlined  in 
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Appendix  A].  The  1/e  width  of  r(Ap)  is  set  to  the  approxi¬ 
mate  width  of  the  Kolmogorov  5/3  power-law  structure  func¬ 
tion,7  the  latter  being  the  more  conventionally  chosen  atmo¬ 
spheric  mode.8  We  take 

r(Ap)  =  exp[-a2(Ap)-].  (2.4) 

where  the  parameter  a  governs  the  width  of  the  wave-front 
correlation  function  relative  to  the  aperture.  This  Gaussian 
correlation  function  was  chosen  for  its  analytical  properties 
and  computational  simplicity.7  A  value  of  a  =  1.3  gives  the 
approximate  match  to  the  Kolmogorov  structure  function. 

Because  the  correlation  function  in  Eq.  (2.4)  is  so  broad 
(when  a  **  1.3),  its  Fourier  transform  is  very  narrow,  and  a 
large  number  of  samples  are  needed  if  one  is  to  avoid  aliasing 
effects.  To  achieve  the  proper  sampling,  we  insert  the  corre¬ 
lation  function  in  a  2046-sample  array;  however,  only  the 
central  128  points  of  each  wave  front  were  used  in  the  correc- 
-tion  tests.  For  the  numerical  computations  carried  out  in 
this  paper,  128  sample  realisations  of  Zip)  were  generated 
(for  fixed  values  of  o 2  and  a). 

Psfs  (<*>)  and  their  corresponding  spatial  transfer  func¬ 
tions 

T(o)  =  |  f(i')e"'“‘de  (2.5) 

were  evaluated,  before  and  after  wave-front  correction,  for 
each  sample  realization  of  the  wave  front  that  was  generated. 

Two  correction  algorithms  were  tested.  In  the  first  ap¬ 
proach,  each  segment  of  the  wave  front  was  fitted  with  a 
straight  line,  using  the  first  and  last  points  of  the  segments 
to  define  the  slopes  of  the  line  segments.  This  is  a  rather 
simplistic  correction  approach,  mimicking  an  algorithm  that 
might  be  used  with  a  segmented  adaptive  mirror.  In  effect, 
it  yields  a  worst-case  result,  since  better  fitting  of  the  seg¬ 
ments  could  be  carried  out  with  more-complex  algorithms. 
The  second  approach  was  to  perform  a  least-squares  polyno¬ 
mial  fit  to  the  wave  front  with  varying  numbers  of  terms  in 
the  fit.  This  approach  would  be  comparable  to  using  a 
rubber  mirror  having  a  limited  number  of  adaptive  degrees 
of  freedom.  We  found  that  when  suitable  parameters  were 
chosen,  the  polynom  al  fitting  gave  similar  results  to  the 
straight-line  fit.  However,  the  Relationship  between  the  fit¬ 
ting  parmeters  and  the  physical  characteristics  cf  the  mirror 
was  far  less  obvious.  Therefore  in  this  paper  we  show  only 
the  results  from  the  straight-line  fit 

3.  RESULTS 

Figure  1  shows  four  of  the  simulated  wave  fronts  generated 
fora  *  1.3  and  (a/ X)  *  2.31,  which  give  a  corresponding  value 
for  (f>/r0)  “  18  85. 7  In  Fig.  2  we  show  examples  of  applying 
straight-line  segment  correction  to  wave  front  A  forO,  2, 4,6, 
and  8  segments.  It  is  clear  that,  even  for  a  few  segments,  the 
phase  excursions  of  the  wave  front  are  dramatically  reduced 
Figu'e  3(a)  plots  the  uncorrected  psf  obtained  from  wave 
from  A.  Figures  3(b),  3(c),  3(d),  3(e),  and  3(f)  show  the  psfs 
obtained  from  wave  front  A  after  straight-line  segment  cor¬ 
rection  has  been  applied  for  2,4,6,  and  8  segment  correction, 
respectively.  Again,  obvious  improvement  in  the  psf  is  ob¬ 
served  even  for  only  two-segment  correction  (the  overall  psf 
width  is  reduced),  and  the  psf  becomes  very  narrow-  with 
increasing  number  of  segments  in  the  correction. 
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In  order  to  test  the  effectiveness  of  the  wave-front  correc¬ 
tion,  we  have  calculated  ensemble-averaged  transfer  fune 
tions  for  the  128  wave-front  realizations  before  and  after 
correction  (Fig.  4),  The  transfer  functions  shown,  after  cor¬ 
rection.  are  for  2,  4,  6,  8,  and  16  straight-line  segment  correc¬ 
tion.  The  16-&egment  correction  result  is  essentially  identi¬ 
cal  to  the  diffraction-limited  transfer  function.  We  also 
have  calculated  the  square  modulus  of  the  ensemble-aver¬ 
aged  transfer  function,  corresponding  to  the  transfer  func¬ 
tion  that  one  would  obtain  for  astronomical  speckle  interfer¬ 
ometry  (Fig.  5).  As  expected,  the  uncorrected  transfer  func¬ 
tion  shows  the  high-frequency  enhancement  expected  in 
speckle  interferometry.  It  is  obvious  from  both  sets  of 
transfer  functions  that  substantial  gains  are  obtained,  even 
when  the  correction-segment  dimension  is  much  longer  than 
the  atmospheric  correlalion  length.  In  most  cases  even  the 
high  frequencies  are  enhanced,  especially  in  the  speckle 
transfer  function. 

The  ensemble-averaged  transfer  functions  were  then  used 
in  a  calculation  of  the  image  for  a  simple  test  object.  An 
object  consisting  of  three  overlapping  Gaussians,  each  hav- 


Apart 


Fig  1.  Four  sample  realizations  of  Zip)  as  measured  in  X  units 
from  an  ensemble  having  a  correlation  function  of  the  form  given  by 
Eq.  (2.4)  with  a  m  1.3  and  <j  •  2.31 


APERTURE 


Fig.  2  Straight-line  corrections  of  order  0.  2,  4, 6,  and  8  applied  to 
sample  realization  A  of  Fig  1 
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Fig  3  Psfs  derived  for  wave-front  realization  A  in  Fig  1.  with 
straight-line  segment  corrections:  lal  no  correction,  (b),  (c),  (dl, 
(e).  (fl  corrections  of  orders  2,  -l.  6. 8.  and  16,  respectively. 


SPATlAv  FREQUENCY 


Fig  4  Square  ensemble  average.  (T(a);\  for  wave-front  realiza¬ 
tion  A  in  Fig  1.  with  straight-line  segment  corrections  of  orders  0  to 

16 

ing  a  different  height  and  width,  was  chosen  to  represent  an 
extended  object  of  interest.  The  Fourier  transform  of  the 
object  was  multiplied  by  various  optical  transfer  functions 
and  then  retransformed  to  image  space.  Figure  6  shows  the 
images  for  the  cases  of  0-,  4-,  8-,  and  16-segment  correction 
along  with  the  diffraction-limited  image 
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The  effectiveness  of  the  correction  was  measured  by  calcu¬ 
lating  the  mean-square  difference  between  the  diffraction- 
limited  image  and  the  uncorrected  and  partially  corrected 
images  In  Figs.  7(a)  and  7(b)  the  results  are  plotted  with 
logarithmic  and  linear  ordinates,  respectively.  We  have 
plotted  both  the  result  in  which  the  average  transfer  func¬ 
tion  was  used  for  the  calculation  and  thi  result  for  wave 
front  C,  one  of  the  least  aberrated  wave  fronts  in  the  ensem¬ 
ble  As  would  be  expected  from  earlier  results,  these  graphs 
show  that  the  corrected  images  are  converging  to  the  unaber¬ 
rated  diffraction  image,  even  for  only  a  limited  number  of 
correction  segments. 

4.  DISCUSSION 

The  numerical  experiments  discussed  in  this  paper  demon¬ 
strate  qualitatively  that  partial  adaptive  correction  of  the 
atmosphere  can  be  an  extremely  useful  technique  lor  astro¬ 
nomical  applications.  However,  the  Gaussian  model  for  the 
atmosphere  is  a  substantial  simplification  and  is  certainly 
not  accurate  at  all  scales  In  addition,  the  effects  of  scintil¬ 
lation,  atmospheric  scattering,  and  telescope  aberrations 
have  not  been  considered.  Experimental  verification  of  the 
results  for  a  range  of  atmospheric  conditions  is  essential 


spatial  frequency 


Fig  5  Speckle  transfer  function.  (lT(o)l,i.  for  wave-front  realiza¬ 
tion  A  in  Fig  1,  with  ttraight-line  segment  corrections  of  orders  0  to 
16. 


Imoge  Plone  Coordmote 

Fig.  6.  Ensemble-averaged  diffraction  images  of  three-peaked 
Gaussian  object  corresponding  to  the  corrected  transfer  functions 
shown  in  Fig.  4 
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Fir.  7.  Integrated  modulus  of  differences  bett  een  the  diffraction 
limited  image  of  the  three-peak  Gaussian  ooject  and  the  corrected 
images  shown  in  Fig  6  (ai  logarithmic  display;  (b)  linear  display. 

before  a  large  effort  is  put  into  implementing  partial  correc¬ 
tion  systems  at  astronomical  sites.  The  range  of  applica¬ 
tions  to  which  the  technique  can  be  usefully  applied  is  also 
highly  dependent  on  the  detailed  atmospheric  parameters. 
Questions  that  may  be  answered  only  by  experimental  mea¬ 
surement  or  a  far  more  detailed  analysis  include  the  follow¬ 
ing:  What  is  the  degree  of  low-frequency  spatial  and  tempo¬ 
ral  correlation?  How  large  is  the  effective  isoplanatic  angle 
at  low  frequencies9  What  is  the  effect  on  the  result  of 
scintillation9  . .  of  scattering9  . .  .  telescope  aberrations9  . .  . 
of  errors  in  the  wn^e-front  sensing  and  correction?  Answer¬ 
ing  these  questions  is  beyond  the  scope  of  the  current  paper, 
requiring  far  more  complicated  modeling  or  experimental 
testing. 

Despite  these  questions,  there  are  many  obvious  astro¬ 
nomical  applications  for  which  partial  adaptive  correction 
could  be  of  substantial  benefit: 

Hi  Many  spectroscopic  sources  are  intrinsically  faint 
and  extended,  requiring  hours  of  integration  time  on  large 
telescopes  Improvements  in  seeing  would  decrease  the  in¬ 
tegration  lime  proportionately,  resulting  in  increased  effi¬ 
ciency  of  large  telescope  time,  a  scarce  and  precious  com¬ 
modity.  For  many  scientific  problems,  the  improved  spatial 
resolution  could  also  be  important.  Of  course,  the  success 
for  faint-object  spectroscopy  will  be  critically  dependent  on 
being  able  to  make  wave  front  measurements  on  off-axis 


reference  stars,  which  csn  be  effective  only  if  the  low-angu¬ 
lar-frequency  isoplantic  angle  is  large. 

(2)  The  search  for  protoplanetary  material  around 
bright  stars  (such  as  the  discovery  of  a  nebula  around  d 
Pictoris9)  requires  coronagraphic  techniques  as  well  as  very' 
low  scattering,  since  the  differences  in  brightness  between 
the  star  and  the  nebula  will  be  large.  Implementation  of 
partial  correction  before  the  coronagraph  not  only  can  per¬ 
mit  more  effective  and  closer-in  blocking  of  the  primary  star 
but  can  also  reduce  the  level  of  scattered  light  in  the  optics. 
Since  the  light  from  the  star  rejected  by  the  coronagraph  can 
be  used  for  wave-front  measurement,  this  should  be  an  ideal 
application  for  correction,  since  no  atmospheric  assump¬ 
tions  are  required  for  its  success. 

(3)  Application  of  speckle  imaging  techiques  after  par¬ 
tial  correction  could  permit  near-real-time  diffraction-limit¬ 
ed  image  reconstruction.  Since  the  number  of  frames  re¬ 
quired  for  convergence  of  the  speckle  process  is  proportional 
to  the  fourth  power  of  the  seeing,  relatively  small  improve¬ 
ments  in  seeing  could  reduce  the  computation  time  required 
for  reconstruction  (without  the  use  of  supercomputers)  suf¬ 
ficient  to  permit  real-time  data  processing  with  convention¬ 
ally  available  hardware. 

It  is  clear  that  partial  correction,  leading  to  improved 
seeing,  can  have  a  significant  effect  on  both  efficiency  and 
ground-based  spatial  resolution  at  many  astronomical  facili¬ 
ties,  since  improved  seeing  is  almost  equivalent  to  an  in¬ 
crease  in  aperture  aize.  Almost  certainly,  as  the  technology 
becomes  available  at  an  affordable  level,  many  new  applica¬ 
tions  will  be  recognized  and  implemented 

APPENDIX  A 

We  briefly  outline  the  procedure  used  to  generate  sample 
realizations  of  the  random  wave  front.  We  let 
,v 

Zip)  m  'S  z„  exp(i0n)exp[i(mT/Lip]  (All 

m  •  -  ;V 

in  the  fundamental  interval  !p  I  <  L.  The  random  variables 
6m  are  uniformly  distributed  in  loj  <  r  and  are  statistically 
independent.  The  z„  are  as  yet  unspecified  coefficients 
Provided  that  N  is  large  enough,  then  Zip)  is  Gaussian 
distributed  by  virtue  of  the  central-limit  theorem. 

We  require  that  the  mean  of  Zip)  vanish: 

N 

(Zip))  -  V  z„(exp(i(?„))exp[i(mir/L)p)  -  0.  (A2) 

m  ■  —  N 

The  only  way  for  this  series  to  vanish  is  for 

<exp(±i0m)  ■  0.  (A3) 

This  is  true  if  the  8„  are  uniformly  distributed  over  (-ir.  *). 

To  evaluate  the  2^,  we  form  the  covariance  of  Zip). 
iZ(Pi)Z*(pj(i;  we  can  show  that 
,v 

(Z(P|)Z*(p3)'  -  V  2„J eip[i (xm/Up].  (A4) 

where  f>  m  \p,  -  p3|.  The  are  given  by 
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=  ~  J*  <Z<p,)Z*(p,)> 

X  exp|i(rrm/Z.)(p,  -  p^jcKp,  -  p2).  (A5) 

For  the  Gaussian  wavs-front  covariance  function,  Eq.  (2.4), 
this  reduces  to 

=  ~  |  expt-.i^lcos^™  pyp.  (A6) 

These  integrals  were  evaluated  numerically.  The  series, 
(Eq  Al),  was  evaluated  by  the  fast  Fourier  transform. 
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ABSTRACT 

High  angular  resolution  observations  of  the  supernova  in  the  Large  Magellanic  Cloud,  SN  1987 A,  have 
revealed  a  bright  source  separated  from  the  SN  by  approximately  60  mas  with  a  magnitude  difference  of  2.7  at 
656  am  (Ha).  Speckle  imaging  techniques  were  applied  to  data  recorded  with  the  CfA  two-dimensional  photon 
counting  detector  (PAPA)  on  the  CTIO  4  m  telescope  on  March  25  and  April  2  to  allow  measurements  in  Ho 
on  both  nights  and  at  533  nm  and  450  am  on  the  second  night  The  nature  of  this  object  is  as  yet  unknown, 
though  it  is  almost  certainly  a  phenomenon  related  to  the  SN. 

Subject  headings  interferometry  —  stars:  superaovae 


1.  INTRODUCTION 

The  supernova  SN  1987 A  in  the  Large  Magellanic  Cloud 
(LMC)  is  the  closest  unobscured  event  of  its  kind  since  1604 
This  provides  a  unique  opportunity  for  a  detailed  study  of 
this  extraordinary  event  using  a  wide  variety  of  observational 
tools.  At  the  50  kpc  distance  of  the  LMC,  high  angular 
resolution  interferometric  techniques  are  capable  of  providing 
substantial  information  about  the  expanding  shell  and  the 
environment  around  the  SN. 

On  1987  March  25  (30  days  after  the  initial  event)  and  on 
April  2,  observauons  of  the  SN  were  made  using  the  Ccrro 
Tololo  Inter  American  Observatory  (CTIO)  4  m  telescope  and 
the  Harvard- Smithsonian  Center  for  Astrophysics  (CfA) 
speckle  imaging  system  The  planned  purpose  of  the  observa¬ 
tions  was  to  attempt  to  measure  the  angular  diameter  of  the 
SN  photosphere  and  to  detect  the  relative  position  of  the  SN 
with  respect  to  the  Sanduleak  B3  supergiant  At  the  lime  of 
these  observations,  it  was  uncertain  whether  this  star  still 
existed,  though  subsequent  analysis  of  IUE  spectral  data  has 
since  shown  (Kirshner  et  a!.  1987)  that  the  star  probably 
cannot  still  be  in  its  original  form.  Speckle  observations  could 
also  be  used  to  examine  the  region  around  the  SN  for  possible 
light  echos  from  surrounding  dust  or  gas.  Despite  its  expected 
angular  size  of  only  a  few  miUiarcseconds  (mas),  measurement 
of  the  photospheric  diameter  appeared  to  be  possible  using 
aperture  mask  interferometry,  since  the  high  brightness  of  the 
source  would  provide  very  high  signal-to-noise  data.  Details 
of  this  experiment  will  be  included  in  a  separate  paper. 

Convenuonal  speckle  data  were  recorded  using  the  PAPA 
two-dimensional  photon  counting  detector  and  a  set  of  nar¬ 
row  band  (10  nm  half-power  width)  interference  filters.  Re¬ 
duction  of  the  data  recorded  at  656  nm  (Ha)  from  both 
nights  produced  a  totally  unexpected  result,  showing  a  bright 
source  only  2.7  mag  fainter  than  the  primary  source,  and 
separated  from  it  by  approximately  60  mas  (Karovska  et  at. 
1987)  A  nearly  identical  result  was  obtained  from  the  Ha 
data  from  the  observations  on  both  nights  and  the  source  was 


also  delected  at  533  nm  and  450  nm  on  the  second  night  This 
second  object  cannot  possibly  be  the  original  Sanduleak  star 
since  it  is  about  6  mag  brighter 

In  this  Letter ,  we  will  describe  in  more  detail  the  observa¬ 
tions  and  results,  suggest  possible  interpretations,  and  discuss 
the  follow-up  observations  needed  to  attempt  to  understand 
the  mysterious  nature  of  this  second  source. 

II.  observations  and  processing 

Two  nights,  1987  March  25  and  April  2,  were  provided  by 
the  director  of  CTIO,  Dr.  R.  Williams,  for  speckle  observa¬ 
tions  of  the  SN.  The  speckle  system  used  for  these  observa¬ 
tions  was  constructed  at  CfA  and  has  been  described  in  detail 
elsewhere  (Karovska,  Nisenson,  and  Siachnik  1986)  An  im¬ 
portant  component  of  this  system  is  the  detector,  the  Preci¬ 
sion  Analog  PhotoD  Address  detector  (PAPA)  (Papaliolios, 
Nisenson,  and  Ebstein  1985)  which  determines  the  x-y  posi¬ 
tion  and  time  of  arrival  of  each  detected  photon  This  detec¬ 
tor  has  been  shown  to  have  the  characteristics  necessary  for 
accurate  speckle  image  reconstruction  The  digital  photon 
addresses  are  converted  to  a  video  signal  and  stored  on  VCR 
tape  for  later  processing.  The  front-end  optics  package  in¬ 
cludes  magnifying  optics  to  match  the  telescope  and  detector 
rescl  scales,  computer-controlled  atmospheric  dispersion  cor¬ 
recting  prisms,  a  set  of  narrow-band  interference  filters,  and 
vanous  neutral  density  filters. 

The  observauons  consisted  of  recording  data  sets  of  5 
minutes  duration  on  the  SN  and  three  different  reference 
stars  Data  were  recorded  with  four  different  fillers  centered 
on  400,  450,  533,  and  656  am,  each  with  a  bandwidth  of  10 
nm  Count  rates  for  the  data  sets  ranged  from  40  to  70 
thousand  detected  photons  per  second  (neuual  density  filters 
were  required  to  reduce  the  source  brightness)  The  field  size 
used  for  recording  was  179  with  256  X  256  pixel  sampling 
This  gives  a  7  mas  per  reset  which  is  more  than  adequate  for 
sampling  the  diffraction  limit  of  the  4  m  telescope,  which 
ranges  from  20  mas  at  400  nm  to  33  mas  at  656  nm  Since  the 
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observations  were  made  when  the  telescope  was  pointed  be¬ 
tween  20°  and  30°  from  the  horizon,  the  data  were  not  of  the 
best  quality.  Seeing  was  approximately  2",  the  atmospheric 
dispersion  was  not  fully  corrected  at  shorter  wavelengths,  and 
the  telescope  tracking  was  unsteady 

The  first  step  in  the  data  processing  is  flat-fielding,  which 
corrects  for  nonuniformities  in  the  photocathodc  and  the 
camera.  The  photons  are  gtouped  into  frames  and  then 
Founer-transformed  (FT)  The  optimum  frame  time  was  de¬ 
termined  by  integrating  short  sets  of  data  with  several  differ¬ 
ent  frame  times  and  measuring  the  vanance  in  the  high 
frequency  tegion  of  the  FT.  The  frame  time  for  the  SN  data 
turned  out  to  be  5  ms  Approximately  60.000  frames  were 
available  for  each  data  set  Complex  correlations  are  calcu¬ 
lated  from  the  FT  of  each  frame  (Knox  and  Thompson  1974) 
to  average  both  the  amplitude  and  phase  in  the  FT;  these 
correlations  are  then  summed  over  all  the  frames.  Similar 
operations  arc  applied  to  the  reference  star  data,  and  division 
of  the  SN  data  by  the  reference  star  data  in  Fourier  space 
performs  an  operation  equivalent  to  deconvolution,  eliminat¬ 
ing  the  effects  of  the  atmosphcnc  transfer  function  and  tele¬ 
scope  aberrations.  The  inverse  FT  of  the  complex  spectrum 
results  in  the  reconstruction  of  an  image  The  transform  of 
the  squared  modulus  reconstructs  an  autocorrelation  image 
(AC);  the  AC  generally  has  better  signal-to-noise  ratio  and  is 
used  for  determining  the  magnitude  difference  of  the  two 
components  In  the  case  of  a  simple  object  such  as  a  binary 
star,  the  image  is  needed  only  to  eliminate  the  180°  ambiguity 
in  the  position  of  the  second  component 

111  RESULTS 

Figure  1  (Plate  LI)  shows  the  results  from  the  data 
processing  Fieure  la  is  a  single  input  frame,  showing  the 
position  of  each  detected  photon  in  a  single  input  frame.  Flat 
fields  were  recorded  by  pointing  the  telescope  at  an  il¬ 
luminated  screen  on  the  dome,  and  the  amplitudes  in  each 
frame  were  scaled  by  the  integrated  flat  field,  thereby  cor¬ 
recting  for  system  and  detector  nonumformities  Figure  16 
displays  the  power  spectrum  from  the  integrated  transform  of 
the  656  nm  data  recorded  on  April  2.  Despite  the  noi:e  in  the 
display,  one  can  see  broad,  low-contrast  bands,  which  are  the 
fnnges  charactenstic  of  a  close  double  star  with  unequal 
magnitude  components  Figure  lc  shows  the  image  of  an 
unresolved  reference  star,  and  Figure  Id  shows  the  SN  recon¬ 
structed  image  The  sue  of  the  displayed  spots  is  indicative 
only  of  their  relative  brightness 

Data  for  nearby  companson  stars  recorded  close  in  time 
through  the  same  filters  produced  clean,  pointlike  images, 
with  no  significant  structure  at  the  separation  and  position 
angle  of  the  second  source  in  the  SN  images.  Figure  2o  plots 
the  Ha  AC,  and  Figure  2 b  shows  the  aC  at  533  am  These 
two  ACs  were  generated  using  the  central  128  x  128  puxels  in 
the  field  (7  mas  sampling,  0795  field)  which  gives  a  somewhat 
better  definition  of  the  reconstructed  peaks.  Figure  2c  shows 
the  AC  at  450  nm,  reconstructed  by  averaging  the  finest  pixels 
to  produce  a  1"9  (128  x  128)  fieid  with  14  mas  sampling 
Figure  2d  plots  the  reference  star  t>  Doradus  which  was 
recorded  with  the  450  nm  data,  also  reconstructed  at  the 
larger  scale.  The  averaging  generally  gives  better  signa)-to- 


noise  ratio,  needed  for  the  450  nm  reconstruction,  but  may 
result  in  an  image  which  is  slightly  less  sharp  It  is  clear  that 
the  double  structure  in  the  SN  is  well  above  the  background 
noise  level,  particularly  in  Ha,  and  the  reference  star  AC  has 
no  substantive  structure  at  a  comparable  scale  and  position. 
The  ACs,  which  are  reconstructed  from  only  the  amplitudes 
in  the  transform,  have  better  signal-to-noise  ratio  than  the 
image,  because  of  the  more  severe  requirements  on  the  quality 
of  the  data  for  the  image  reconstruction  process  Therefore, 
they  were  used  to  measure  the  parameters  of  the  recon¬ 
structed  sources. 

The  separation  of  the  two  sources  measured  in  Ha  on  both 
nights  was  0"059  ±  0"008.  and  the  magnitude  difference  was 
2.7  +  0.2.  Since  the  visual  magnitude  of  the  SN  at  the  time  of 
the  observations  ranged  from  4  0  (March  25)  to  3.8  (April  2). 
this  source  has  an  apparent  magnitude  of  about  6  5  The 
errors  in  our  magnitude  estimates  are  slightly  too  large  (0  2) 
to  allow  a  determination  whether  the  second  source  changed 
in  brightness  by  the  same  amount  as  the  SN  during  the  eight 
days  between  our  observations  The  measurement  at  533  nm 
(April  2)  gave  a  separation  of  0"052  ±  0'007,  and  the  mea¬ 
sured  magnitude  difference  was  3.0  ±  0.5.  In  both  the  Ha 
and  533  nm  measurements,  the  position  angle  of  the  sec¬ 
ondary  source  relative  to  the  primary  was  194°  +  2°  The 
180°  ambiguity  normally  associated  with  speckle  interferome¬ 
try  measurements  was  eliminated  using  the  image  reconstruc¬ 
tion  Reconstructions  from  data  recorded  at  450  nm  show  a 
feature  at  approximately  the  same  position  with  a  3. 5-4.0 
mag  difference  from  the  primary  Residual  atmospheric  d's- 
persion  produced  elongation  in  the  reconstructions,  and  this 
was  an  increasingly  severe  problem  at  shorter  wavelengths 
Thus  either  residua)  dispersion  or  the  red  color  of  the  source 
might  explain  the  uncertainty  of  its  detection  at  450  nm  The 
elongation  appears  to  be  somewhat  greater  in  the  SN  images 
than  in  the  comparable  reference  stars,  despite  their  having 
been  close  in  position  or,  the  sky  This  suggests  that  the 
elongation  may  be  in  the  object,  though  the  signal-to-noise 
ratio  in  the  detection  is  loo  low  to  be  sure 


tv.  DISCUSSION 

The  detection  of  the  second  bright  source  on  two  different 
nights.  8  days  apart,  combined  with  the  reported  detection  of 
this  source  by  a  group  at  Imperial  College  performing  speckle 
observations  on  the  Anglo-Auslralian  Telescope  (Matcher, 
Meikle,  and  Morgan  1987)  on  April  14,  leaves  little  doubt 
that  the  source  really  exists  Since  this  second  source  was 
almost  5  mag  brighter  than  any  known  preexisting  source  in 
the  field  (the  B3  Sanduleak  star),  it  r  clear  that  its  appearance 
must  be  related  to  the  SN.  Reconstructions  from  the  data 
recorded  at  400  nm,  where  the  SN  was  reduced  sufficiently  in 
brightness  so  that  the  magnitude  difference  between  the  SN 
and  the  B3  star  would  have  been  small  enough  to  delect  their 
separation,  showed  no  extension  of  the  primary  source  to  an 
accuracy  of  about  20  mas.  The  60  mas  separation  of  the  two 
bright  sources  would  correspond  to  about  3000  AU  at  the 
distance  of  the  LMC  (50  kpc),  perpendicular  to  the  line  of 
sight  This  corresponds  to  2  it-weeks,  and  the  first  observation 
occurred  30  days  after  the  SN  explosion 
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One  possible  explanation  for  the  existence  of  the  second 
source  is  that  it  could  be  a  large  cloud  of  dust  or  gas  reflecting 
light  from  the  SN.  However,  the  brightness  of  the  second 
source  would  require  that  the  cloud  subtend  a  minimum  of 
1  sr  as  seen  from  the  SN  (the  source  is  10%  as  bright  as  the 
SN),  even  with  an  albedo  of  unity  This  would  be  an  enor¬ 
mous  cloud,  probably  spatially  resolved  by  our  measure¬ 
ments  It  is  also  unlikely  that  such  a  cloud  could  exist  in  the 
neighborhood  of  a  B3  supergiant  Another  possibility  is  that 
the  second  source  was  a  very  bnght  star  wrapped  in  a  dust 
shell  which  was  blown  away  by  the  SN  flash.  However,  such  a 
source  would  be  nearly  the  brightest  star  in  the  LMC  (im¬ 
probable)  and  would  have  been  very  bnght  in  the  IR  No  such 
source  existed  in  that  position  in  the  IRAS  survey  of  infrared 
sources 

While  there  are  many  other  more  exotic  possibilities  to 
explain  these  results,  it  is  clear  that  more  data  are  needed  to 
sort  them  out  Measurement  at  a  later  date  should  determine 
whether  their  relative  positions  have  remained  stationary  and 
whether  the  second  source  has  brightened  along  with  the  SN 


Measurement  at  other  wavelengths  should  determine  what  the 
relative  color  of  the  two  sources  is,  and  whether  the  second 
source  is  an  emission-line  object  New'  observations  are 
planned  at  the  end  of  1987  May  and  in  early  July.  It  is  hoped 
that  they  will  shed  new  light  on  this  extraordinary  object. 
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4  m  observing  time  on  such  short  notice,  and  to  Oscar  Saa 
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also  wish  to  thank  N.  Carleton  and  S.  Ebstein  for  their  aid  in 
our  frantic  last-second  preparations  of  the  equipment;  1  1 
Shapiro  for  his  rapid  support  of  the  expedition;  and  R 
Kjrshner,  J  Raymond,  L  Hartmann,  and  W.  Traub  for 
numerous  useful  discussions  and  advice  We  are  also  grateful 
to  Henry  Radoski  of  AFOSR  for  his  continued  belief  in  and 
support  of  the  speckle  program  at  CfA.  This  work  has  been 
supported  under  grant  AFOSR-86-0103.  NASA  grant  NGL- 
22-007-228,  and  Smithsonian  Institution  Research  Opportuni¬ 
ties  and  Scholarly  Studies  Programs 
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THE  supernova  SN1987A  io  the  Large  Magellaalc  Cloud  has  been 
observed  by  high-angular-resolutioo  speckle  interferometry  since 
25  March  (30  days  after  the  explosloo)  with  the  4-m  telescope  at 
the  Cerro  Tololo  loteramerican  Observatory.  These  observations 
have  provided  a  number  of  results  that  may  be  central  to  a  detailed 
understanding  of  this  unique  event.  Data  obtained  on  25  March 
and  2  April  1987  revealed  a  second  bright  'companion'  source 
separated  from  the  rupernova  by  60  tnilliarcsecoods  and  less  than 
three  magnitudes  fainter  than  the  supernova1.  Measurements  of 
the  average  diameter  of  the  supernova  envelope  have  been  made 
from  data  recorded  from  March  1987  to  April  19882.  Here  we 
preseot  a  more  detailed  analysis  of  these  data,  which  ahows  that 
the  expanding  envelope  is  asymmetric.  This  asymmetry  has  been 
observable  since  June  15*87.  The  ratio  between  the  minor  and  major 
axes  of  the  envelope  profile  Is  about  2-3,  and  the  position  angle 
of  the  major  axis  is  20°  ±  5°,  consistent  with  results  reported  from 
polarization  measurements.  The  major  axis  is  aligned  with  the 
position  angle  of  the  companion  to  the  supernova. 

Speckle  observations  at  optical  wavelengths  have  been  earned 
out  using  the  PAPA  detector’,  a  two-dimensional  photon- 
counting  sensor  which  records  a  catalogue  of  successive  photon 
positions  The  detector  currently  has  a  maximum  data  recording 
rate  of  100,000  photons'1  in  a  field  of  512x512  pixels.  The 
speckle  camera  uses  a  fore-optics  package  which  provides  mag¬ 
nification  of  the  image,  so  that  the  diffraction-limited  scale  of 
the  telescope  is  matched  to  the  pixel  size  of  the  camera,  narrow- 
band  optical  filtering,  which  yields  temporal  coherence  sufficient 
to  allow  interference  over  the  path-length  errors  introduced  by 
atmospheric  aberrations,  and  atmospheric-dispersion  correc¬ 
tion,  using  a  computer-controlled  prism  compensator.  The 
digital  output  of  the  camera,  in  the  form  of  photon  addresses, 
is  encoded  on  a  video  earner  and  recorded  with  a  VCR.  The 


stored  addresses  are  recovered  in  the  laboratory  and  unpacked 
for  the  speckle  processing. 

The  photon  addresses  are  grouped  into  frames,  with  a  short 
time  per  frame  relative  to  the  atmospheric  correlation  time.  This 
approach  allows  framing  of  the  data  for  maximization  of  the 
stgnal-to-noise  ratio  in  the  integrated  result;  typical  frame  times 
range  from  2  to  20  ms.  individual  frames  are  built  from  the 
photon  list  by  incrementing  the  number  in  the  array  position 
corresponding  to  the  address  of  each  detected  photon.  Correc¬ 
tions  for  the  camera  sensitivity  (flat  fielding)  are  made  on  the 
resulting  frames.  The  Fourier  transform  for  each  frame  is  then 
calculated  and  accumulated  into  the  power  spectra  and  the 
complex  correlation  arrays  which  are  required  for  diameter 
measurement  and  speckle  image  reconstruction  using  our  ver¬ 
sion  of  the  Kjiox-Thompson  image-reconstruction  algorithms4. 
Compensation  for  the  atmospheric  and  telescope  transfer  func¬ 
tion  is  accomplished  by  observing  an  unresolved  comparison 
star  immediately  before  and  after  the  object  observations  This 
star  is  chosen  to  be  as  close  in  angular  position  to  the  object  as 
possible  so  that  its  atmospheric  statistics  are  similar.  Deconvol¬ 
ution  by  the  reference  star  enhances  the  amplitudes  of  the  high 
angular  frequencies  in  the  reconstruction.  A  detailed  treatment 
of  the  algorithms  and  data  processing  is  discussed  elsewhere4. 

Data  sets  on  SN1987A  and  other  stars  result  in  diameter 
determinations  with  a  precision  of  better  than  one  milliarcsec. 
These  measurements  are  made  by  fitting  the  integrated  power 
spectra  to  the  power  spectrum  of  a  uniform  disk1,  and,  depend¬ 
ing  on  the  signal-to-noise  ratio,  can  give  results  with  resolution 
exceeding  the  diffraction  limit  of  the  telescope.  Further  analysis 
of  the  speckle  data  recorded  between  June  1983  (95  days  after 
the  explosion)  and  Apnl  1988  (411  days  after  the  explosion) 
shows  that  the  expanding  shell  is  elongated5  This  asymmet' 
is  detected  from  data  recorded  ai  several  wavelengths  between 
442.0 nm  and  850.0  nm.  Figures  la,  c,  t  and  g  show  the  power 
spectra  of  SN1987A  obtained  from  the  data  recorded  near  the 
centre  of  the  Ha  emission  line  in  May-June  >987,  November 
1987,  February-  March  1988  and  April  1988.  These  power  spectra 
are  elongated,  with  a  major  axis  corresponding  to  a  position 
angle  (PA)  in  image  space  of  20°  (or  200°)i5°.  Similar  depar¬ 
tures  from  circular  symmetry  are  seen  in  the  power  spectra 
obtained  from  data  recorded  in  several  different  wavelengths. 
By  fitting  the  power  spectra  of  uniform-bnghtness  ellipses  to 
the  observed  power  spectra,  we  determined  the  lengths  of  the 


FIG  1  Speckle  power  spectra  ot  SN1987A  recor¬ 
ded  m  Ho  at  four  Oiterent  epochs  (fl  c  e  g)  and 
ellipses  snowing  schematically  the  sue  and 
orientaiion  of  the  enveiooe  at  each  epoch  10.  A  f 
h)  a.  0.  June  1987.  c.  o.  November  1987:  e.  f. 
February  1988  &.  h  April  1988 


major  and  minor  axes,  and  their  position  angles.  All  or  the 
measurements  gave  the  same  PA  within  the  error  bars.  A  sum¬ 
mary  of  the  measurements  is  given  in  Table  1  In  contrast,  data 
recorded  on  comparison  stars  give  symmetric  power  spectra  and 
symmetric  images  with  measured  angular  diameters  of  less  than 
3  mas. 

The  ratio  of  the  minor  to  major  axis  of  the  supernova's 
envelope  is  about  2-3  for  all  four  data  sets.  The  apparently 
increasing  asymmetry  in  the  supernova's  power  spectra  shown 
in  Fig.  1  is  due  to  the  increasing  angular  size  of  the  supernova 
relative  to  the  telescope’s  diffraction  limit.  Figures  1  i>,  d,  /  and 
h  show  the  ellipses  whose  Fourier  transforms  are  fitted  to  the 
power  spectra  at  the  vanous  epochs. 

True  images  of  SN1987A  and  its  comparison  star,  v  Vol.  were 
reconstrt'cted  using  Knox-Thompson  speckle-imaging  tech¬ 
niques'*.  Figure  2  shows  the  images  from  data  recorded  in  April 
1988  at  533.0  nm  (IO-nm  bandpass)  and  Ha  (10-nm  bandpass). 
The  elongation  of  the  supernova  images,  along  an  axis  inclined 
at  a  position  angle  of  20°  (or  200°),  is  quite  evident  despite  the 
fact  that  its  angular  size  in  April  was  only  slightly  greater  than 
the  telescope's  diffraction  limit.  Because  the  size  of  the  image 
was  so  close  to  the  telescope  diffraction  limit,  the  size  of  the 
axes  were  determined  from  power  spectra,  not  from  recon¬ 
structed  images.  The  images  also  appear  to  be  somewhat  brighter 
in  the  south-west  direction.  The  x-  Vol  images  show  no  elongation 
and  are  symmetric.  Their  size  is  equivalent  to  the  beam  size  for 
the  reconstruction  process. 

In  addition  to  our  speckle  results,  other  observations  of 
SN1987A  show  strong  evidence  for  an  asymmetrically  expand¬ 
ing  shell*’.  Early  spectroscopic  observations  (between  20  and  80 
days  after  the  explosion)  detected  a  double-peaked  feature  in 
several  hydrogen  lines  at  optical  and  infrared  wavelengths7-1. 
The  appearance  of  two  equally  displaced  components  at  longer 
and  shorter  wavelengths  than  the  maximum  emission  was  inter¬ 
preted  as  being  due  to  some  departure  from  spherical  symmetry 
of  the  shell.10  Polarimetric  observations1114  also  show  strong 
evidence  for  asymmetry  in  the  envelope  with  a  position  angle 
of  -200°.  The  observed  polarization  structure  was  interpreted 
as  arising  from  an  asymmetrically  expanding,  scattering  atmos¬ 
phere.  Polanmetnc  data  were  modelled  by  a  scattering  almos- 


riC  2  Image  reconstructions  of  SN1987A  and  a  comparison  star  recorded 
in  April  1988  a  SN1967A  at  533  0  nm  c.  Comparison  star  at  533  0  nm 
c  SN1987A  at  ha.  a  Comparison  star  at  ha 


TABLE  1 

Measured  asymmetries  for  SN1987A  (PA 

-200  *  5°) 

Days  after 

Average 

diameter 

Major 

axis 

Minor 

axis 

explosion 

A  (nm) 

(mas) 

(mas) 

(mss) 

95-98 

5330 

18  x  2 

19*2 

17*2 

(May-June  87) 

6  565  (Ha) 

3*1 

10*2 

6*2 

7,750 

15*2 

20*2 

10*2 

265-268 

6.585  (Ha) 

17*1 

18*3 

14*3 

(Nov.  87) 

8500 

16*1 

20*2 

14*2 

370-373 

4420 

25*2 

27*2 

22*2 

(Pe6  -Mar  88) 

5330 

20*1 

24*2 

14*2 

6400 

17*1 

20*2 

15*2 

6.565  (Ha) 

21*1 

24*2 

16*2 

7,000 

17*1 

20*2 

16*2 

8500 

18*2 

25*2 

12*2 

409-411 

5330 

26*2 

30*2 

20*2 

(April  88) 

656S(Ha) 

27  il 

30*2 

20*2 

8.500 

24  ±2 

28*3 

16*3 

phere  shaped  as  a  prolate  or  oblate  spheroid,J1}'l\  Best  fits 
have  been  obtained  for  ratios  of  the  smallest  to  the  largest  axis 
of  these  spheroids  ranging  from  0.6  to  0.9.  Recently,  y-ray 
observations  have  yielded  linewidths  requiring  either  fragmenta¬ 
tion  or  asymmetry  of  the  supernova  shell  for  their  explanation10. 

Fhesent  results  from  speckle  interferometry,  polarimetry,  spec¬ 
troscopy  and  y-ray  observations  show  that  the  expanding  atmos¬ 
phere  of  SN1987A  is  not  spherically  symmetric.  The  theoretical 
models  that  assume  spherical  symmetry  must  clearly  be  revised. 
An  asymmetric  envelope  may  also  be  necessary  to  explain,  for 
example,  the  early  emergence  of  X-rays  and  y-rays  and  the 
evolution  of  their  spectra  with  time17.  Woosley1*  has  suggested 
that  non-uniformities  in  the  collapse  and  core  bounce  could 
produce  an  asymmetric  explosion,  but  that  this  asymmetry  is 
far  less  likely  to  propagate  to  the  outer  envelope.  Chevalier  and 
Soker”  have  modelled  the  expanding  supernova  envelope  and 
concluded  that  the  most  likely  mechanism  for  asymmetry  is 
rotational  flattening  of  the  progenitor  envelope,  probably  requir¬ 
ing  a  binary  companion  during  stellar  evolution.  Much  of  the 
evidence  for  asymmetries  and  extended  structures  around  the 
supernova  is  summarized  by  Trimble20.  □ 
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ABSTRACT 

We  present  direct  measurements  of  the  angular  diameter  of  SN  1987A  in  the  Large  Magellanic  Cloud 
(LMC)  made  from  high-angular  resolution  observations  at  the  CTIO  4  m  telescope  in  five  observing  runs 
from  1987  April  to  1988  April.  Diameters  were  determined  to  milliarcsecond  precision  from  integrated  power 
spectra,  using  speckle  interferometric  data.  The  accuracy  of  the  technique  was  evaluated  by  laboratory 
experiments  and  measurements  of  the  diameters  of  several  stars  of  known  angular  size.  The  SN  1987A 
diameters  measured  near  the  center  of  the  Ha  line  show  clear  evolution  during  this  period.  The  rate  at  which 
the  supernova  size  changed  at  this  wavelength  corresponds  to  2850  km  s'1  mean  velocity  of  expansion. 
Diameter  measurements  obtained  in  several  spectral  lines  and  in  the  continuum  indicate  stratification  of  the 
expanding  envelope  of  the  supernova.  Our  continuum  data  yield  diameters  substantially  larger  than  those 
calculated  from  photometric  measurements  and  a  blackbody  fit  to  the  observed  spectra 
Subject  headings:  interferometry  —  stars:  individual  (SN 


I.  INTRODUCTION 

High-angular  resolution  observations  of  SN  1987A  provide 
valuable  data  that  will  contribute  to  the  understanding  of  this 
extraordinary  astrophvsical  event.  Using  speckle  interferomet¬ 
ric  techniques,  it  is  possible  to  overcome  the  resolution  limits 
imposed  by  atmospheric  turbulence  and  observe,  in  some 
detail,  the  morphology  of  the  supernova  (SN)  and  its 
environment.  A  wide  variety  of  other  observational  techniques 
such  as  spectroscopy,  photometry,  and  polanmetry  are  provid¬ 
ing  an  extensive  data  base  that,  combined  with  the  high- 
angular  resolution  measurements,  will  result  in  a  detailed 
understanding  of  the  mechanisms  and  dynamics  of  super- 
novae,  thereby  testing  the  validity  of  theoretical  models. 

Speckle  observations  at  CTIO  in  1987  March  and  April 
(Nisenson  et  al.  1987)  have  already  revealed  the  presence  of  a 
second  bright  source  60  milliarcseconds  (mas)  from  the  SN  and 
within  a  factor  of  12  of  its  brightness.  Observations  at  the  A  AT 
in  1987  April  (Meikle,  Matcher,  and  Morgan  1987),  also  using 
speckle  techniques,  produced  supporting  evidence  for  the 
second  source.  Furthermore,  analysis  of  data  acquired  in  five 
observing  runs  between  1987  April  and  1988  April  have  pro¬ 
vided  measurements  of  the  apparent  diameter  of  the  SN  at 
several  wavelengths  Diameter  measurements  were  made  by 
fitting  the  integrated  power  spectra  obtained  from  the  speckle 
process  to  the  power  spectrum  of  a  uniform  disk  or  a  limb- 
darkened  disk.  Accurate  fitting  allows  diameter  estimates  at 
scales  well  below  the  “diffraction  limit"  of  the  telescope,  to 
a  precision  of  a  few  milliarcseconds  Figure  I  illustrates  the 
principle  by  which  diameters  well  under  the  diffraction  limit  of 
the  imaging  system  (//£>,  where  /.  is  the  wavelength  and  D  is  the 
telescope  diameter)  can  be  measured,  showing  the  power 
spectra  for  stars  with  diameters  0.5,  1,  and  2  times  the  diffrac¬ 
tion  limit.  A  star  having  an  angular  diameter  smaller  than  the 
diffraction  limit  has  a  partial  drop  in  power  at  the  highest 
frequencies  measured.  Accurate  measurement  of  this  drop 
allows  an  estimation  of  the  stellar  diameter  and  this  accuracy  is 
only  limited  by  the  signal-to-noise  ratio  in  the  data.  A  recent 

1  Visaing  Sciemisi  from  CERGA,  Avenue  Copernic.  061 30Grssse,  France 
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demonstration  of  this  principle  has  been  given  by  Davis  and 
Tango  (1986).  Using  a  two-element  interferometer  with  an  1 1 
m  baseline  which  has  a  resolution  limit  of  10  mas  (at  550  nm). 
they  were  able  to  measure  the  diameter  of  Sirius  (5.5  mas)  to  a 
precision  of  0.08  mas.  Similar  results  have  been  obtained  using 
the  I2T  interferometer  at  CERGA  (Koechlin  and  Rabbia 
1985).  While  measuicments  of  equivalent  precision  have  not 
previously  been  demonstrated  using  speckle  techniques,  the 
validity  of  our  results  is  supported  by  laboratory  simulations 
and  measurements  of  stars  whose  diameters  have  been  deter¬ 
mined  by  other  methods. 

it.  observations  and  process.ng 

We  report  results  from  five  observing  runs  using  the  Cerro 
Tololo  Inter- American  Observatory  (CTIO)  4  m  telescope.  A 
summary  of  the  observations  is  given  in  Table  1  Data  were 
recorded  using  the  PAPA  two-dimensional  photon  counting 
detector  (Papaliolios,  Nisenson,  and  Ebstein  1986)  and  a  front- 
end  optics  package  The  optics  package  includes  magnifying 
optics  to  match  the  telescope  and  detector  pixel  scales,  atmo¬ 
spheric  dispersion-correcting  pnsms,  a  set  of  narrow-band 
interference  fillers,  and  various  neutral  density  filters.  The 
camera  records  the  x-y  position  and  time  of  arrival  of  each 
photon  event  as  it  is  detected,  allowing  construction  of  short- 
exposure  speckle  frames  during  the  computer  processing  oper¬ 
ations.  A  major  advantage  of  this  approach  is  that  the 
optimum  exposure  time  (as  determined  by  the  atmospheric 
correlation  time.)  that  maximizes  the  signal-to-noise  ratio  in 
the  reconstruction  may  be  determined  dunng  the  data  pro¬ 
cessing. 

Data  were  recorded  at  many  different  wavelengths  and 
under  a  variety  of  atmospheric  conditions.  The  general 
approach  was  to  record  a  relatively  short  data  set  (5  or  10 
minutes)  on  the  SN,  preceded  and  followed  by  a  5  minute  set 
on  a  comparison  star.  In  all  cases,  the  comparison  stars  were 
chosen  to  be  single  stars  close  in  angular  position  to  the  SN 
and  to  have  expected  angular  sizes  as  small  as  possible  (less 
than  1  mas).  The  comparison  stars  were  also  chosen  to  be  as 
bright  as  possible,  yielding  maximum  signal-to-noise  ratio  in 
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Fic  1  —Ideal  power  specira  for  stars  with  diameters  (<sl  0  5.  (6 I  I.  and  (cl  2  times  the  telescope  diffraction  limit  The  right-hand  axis  corresponds  to  the  4  m 
diffraction  limn  at  6000A  (33  arcsec"  ‘] 


the  time  alloted  for  their  observation,  since  they  are  used  in  a 
deconvolution  calculation  during  the  data  processing-  The  dia¬ 
meter  of  the  detector  field  for  all  data  sets  was  1  '.'95  with 
256  x  256  sampling,  giving  a  7.5  mas  pixel  size.  The  SN  and  its 
comparison  star  were  never  higher  than  an  altitude  of  55°,  and 
were  as  low  as  20°  for  some  of  our  observations,  so  they  all 
required  a  large  and  accurate  correction  for  atmospheric  dis¬ 
persion  Precise  correction  was  performed  using  the  pairs  of 
computer-controlled  Risley  prisms  in  our  optical  front  end 

TABLE  I 


Summary  or  Speckle  Observations 


Date 

Days  after 
Explosion 

Central 

Wavelength 

(nm) 

Bandpass 

(FWHM) 

(nm) 

198?  Apr  2  . 

38 

4502 

10.0 

532.9 

8.0 

656.5 

107 

198?  May  30-Jun  2  . 

95-98 

450  0 

100 

5329 

8.0 

6565 

107 

7750 

100 

1987  Nov  15-18 

265-268 

532  9 

8.0 

640.0 

10.0 

656  5 

107 

658  5 

too 

700.0 

100 

850  0 

250 

1988  Feb  29-Mar  3  . 

370-373 

4420 

32  8 

5330 

80 

6400 

10.0 

656.5 

107 

700.0 

100 

850.0 

250 

1983  Apr  8-10  . 

409-411 

5330 

80 

656  5 

107 

850  0 

250 

Since  the  comparison  star’s  dispersion  will  be  almost  identical 
to  that  of  the  SN,  any  residual  error  due  to  dispersion  will  be 
corrected  in  the  deconvolution  during  the  data  processing.  The 
fit  was  only  performed  out  to  frequencies  where  both  the  refer¬ 
ence  and  object  power  spectra  were  nonzero.  The  seeing  was 
generally  2"-3"  except  in  1988  February.  March,  and  April 
when  seeing  approached  1". 

The  first  step  in  data  processing  is  to  calculate  a  flat  field 
that  corrects  for  the  nonuniformities  present  in  the  photo- 
cathodc  sensitivity  and  the  camera  optics.  The  photons  are 
grouped  into  speckle  frames  and  each  frame  is  flat-fielded 
before  being  Fourier  transformed.  The  procedure  for  flat- 
fielding  the  data  (Ebstein  1987)  is  to  first  generate  a  flat-field 
"mask”  using  a  combination  of  “ white-field ”  data  (in  which 
the  camera  is  looking  at  the  dome  or  an  internal  uniform 
source)  and  a  filtered  version  of  the  long-exposure  image  of  the 
object  being  flat-fielded  Use  of  the  object  data,  itself,  allows  for 
any  dynamic  changes  in  the  flat  field  and  reduces  the  require¬ 
ment  of  recording  “white-field"  data  to  about  once  per  night 
The  synthesized  flat  field  is  then  applied  to  each  frame,  after  the 
photon  addresses  have  been  binned  into  individual  frames. 
One  ends  up  with  frames  in  which  each  photon  is  given  a 
fractional  height  proportional  to  the  flat-field  amplitude  at 
that  position. 

The  optimum  frame  time  is  determined  by  integrating  the 
comparison  star  power  spectra  for  several  different  frame 
times,  calculating  the  variance  and  signal-to-noisc  ratio  in  the 
high-frequency  region  of  the  power  spectrum  and  choosing  the 
exposure  time  which  maximizes  the  signal-to-noise  ratio. 
Typical  frame  times  were  found  to  be  about  10  ms,  though  they 
varied  between  5  and  1 5  ms. 

The  ensemble  average  power  spectrum  is  calculated  for  both 
the  SN  and  its  comparison  star.  Division  of  the  SN  power 
spectrum  by  the  comparison  star  power  spectrum  corrects  for 
the  effects  of  the  atmospheric  and  the  telescope  transfer  func¬ 
tions.  Diameter  measurements  were  made  by  first  azimuthally 
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averaging  the  power  spectra  to  enhance  the  signal-to-noise 
ratio  before  division.  A  visibility  function  (the  power  spectrum 
of  a  limb-darkened  or  a  uniform  disk),  in  the  form  of  the 
Hankel  transform  of  the  stellar  profile  was  least-square-fittted 
to  the  ratioed  power  spectra  with  two  free  parameters,  the 
width,  a,  and  the  amplitude,  m.  The  diameter  in  image  space 
was  calculated  from  the  width  parameter,  a.  In  the  fitting  pro¬ 
cedure,  the  data  were  weighted  by  the  autocorrelation  of  a 
circle  This  corresponds  to  the  fall  off  in  signal  to-noise  ratio 
with  increasing  spatial  frequency  of  an  ideal  telescope  with  no 
atmosphere.  Errors  were  estimated  by  calculating  the  standard 
deviations  of  the  data  from  the  fits.  When  the  supernova  power 
spectrum  was  divided  by  the  corresponding  comparison  star 
power  spectrum,  the  signal-to-notse  ratio  at  the  higher  fre¬ 
quencies  in  both  object  and  comparison  star  power  spectra 
was  low,  so  fitting  was  restricted  to  the  low  and  intermediate 
frequency  range.  In  addition,  fluctuations  in  seeing  between 
recording  of  the  object  and  the  comparison  star  make  correc¬ 
tion  at  the  lowest  few  frequencies  inaccurate.  Despite  these 
restrictions,  this  approach  appears  to  give  reproducible  results, 
and  application  of  this  technique  to  stars  with  known  diameter 
gives  accurate  results. 

in.  laboratory  simulations 

To  check  this  technique,  the  same  modeling  and  fitting  was 
applied  to  data  generated  in  laboratory  simulations  with  disks 
of  known  diameter.  The  speckle  camera  system  (including  the 
PAPA  detector)  was  set  up  to  accept  light  from  an  optical 
system  which  simulates  a  telescope  and  a  randomly  varying 
atmosphere.  The  laboratory  optics  were  scaled  so  that  a  100 
pinhole  p  i  85  cm  from  the  telescope  aperture  was  just 
resolved  at  65o  urn.  The  diameter  of  the  optical  system  pupil 
was  adjusted  so  that  the  moving  piece  of  ground  glass  used  to 
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Fig  2  —  Laboratory  simulation  power  spectra  and  beit-fil  Airy  daks  for  the  50  and  100  jjm  pinholes  Vertical  dolled  lines  show  the  range  of  filled  region  The 
diffraction  lirail  of  the  opucaJ  system  at  650)  A  was  30  arcaec " 1 . 


simulate  atmospheric  turbulence  produced  speckle  images 
with  the  equivalent  of  l"-2”  seeing.  The  light  level  and  spectral 
passbands  were  set  to  match  the  SN  observation  levels. 
Speckle  data  were  recorded  for  three  pinholes  of  diameters 
12  +  2,  50  ±  5,  and  100  +  5  /im;  the  12  jtm  pinhole  was  used 
as  the  comparison  source  after  correction  for  its  finite  size  by 
division  of  the  appropriate  Airy  function.  The  power  spectra 
and  the  best-fit  Airy  functions  for  the  50  and  100  pm  data  are 
shown  in  Figure  2.  The  diameters  calculated  from  the  best  fits 
to  the  data  were  55  +  3  and  100  +  3  gm  The  accuracy  of  these 
results  strongly  supports  the  validity  and  precision  of  this  tech¬ 
nique  for  measuring  diameters  smaller  than  the  diffraction 
limit  of  the  imaging  optical  system. 

iv.  results 

In  an  important  test  of  the  performance  of  our  technique,  we 
measured  the  diameter  of  a  solar-type  star,  a  Cen  A.  The  data 
were  recorded  using  the  CTIO  4  m  telescope  on  1987  June  1  ai 
533.0  nm  and  1988  February  29  at  450.0  nm.  From  both 
observing  runs  we  obtained  9  +  2  mas  for  the  angular  dia¬ 
meter,  assuming  a  uniform  disk  Figure  3  shows  the  data  from 
1988  February  and  the  corresponding  best  fit  When  the  effect 
cf  limb  darkening  is  taken  in  account  (Allen  1973,  p.  120)  the 
measured  diameter  of  a  Cen  A  increases  by  approximately 
10%.  Our  diameter  measurement  for  a  Cen  A  is  very  close  to 
the  diameter  of  8.6  +  0.2  mas  obtained  photometrically  by 
Blackwell  and  Shallis  (1977)  The  diameter  of  the  red  super¬ 
giant  a  Sco  at  533.0  and  656.5  nm  was  also  measured.  Figure  3 
shows  the  data  recorded  in  1988  April  (at  533.0  nm)  using  the 
CTIO  4  m  telescope  and  the  corresponding  best  fit.  At  the  two 
wavelengths  of  observation  we  obtained  35  +  1  mas  for  the 
diameter  of  a  Sco,  assuming  uniform  disk.  When  limb  dark¬ 
ening  is  included,  assuming  a  linear  limb-darkening  law,  a  dia- 
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Fig  3  —Power  specira  for  i  Cen  A  and  a  Sco.  surs  with  known  diameters.  The  best  fits  to  the  data  (corresponding  to  a  uniform  disk  of  9  mas  for  a  Cen  A  and  33 
mas  for  a  Sco)  are  also  plotted  Vertical  dotted  lines  show  the  range  of  fitted  region 


meter  of  39  ±  1  mas  was  obtained.  These  measurements  are 
very  close  to  the  diameter  of  36  ±  3  mas,  obtained  by  Black- 
well  and  Shallts  (1977).  Also,  they  are  in  good  agreement  with 
the  measurements  from  several  observers  using  different  tech¬ 
niques  (summarized  by  White  1980) 

Table  2  shows  the  summary  of  the  results  of  diameter  mea¬ 
surements  of  SN  1987A  (assuming  a  uniform  disk)  made  from 
data  recorded  at  the  five  different  epochs  corresponding  to  the 
38th.  95th-98th,  265th-268th,  370th-373d.  and  409th— 41  Ith 
day  after  the  explosion  While  the  data  are  not  of  sufficient 
quality  to  allow  us  to  calculate  the  profile  of  the  SN,  the  effect 
of  iimb  darkening  (or  brightening)  on  the  measured  diameters 
can  be  estimated  Calculations  show  that  extreme  limb  dark¬ 
ening  (or  brightening)  will  change  the  diameter  measurement 
by  no  more  than  15%. 

From  the  March-April  observations  (at  wavelengths  of 
450.0,  533.0.  and  656.5  nm).  the  SN  was  found  to  be  resolved  at 


TABLE  2 

Angular  Diameter  Measurements  of  SN  I987A* 


Days  after  Explosion 

A  - - - — — — 

mi)  38  9S-98  265-268  370-373  409-411 


442  0  ...  25  +  2 

4500  12  T  5  23  :  4 

532  9  11+  4  18  ;  2  20  +  2  20  ±  1  26  +  2 

640  0  . ....  15  ;  1  17+1 

656  5  S5  8  ;  I  18+1  21  ±  I  27+1 

658  5 .  17  ±  I 

7000 .  17  +  I  17+1 

775  0 .  ...  15  +  1 

850.0 .  16 ;  I  18  +  2  24  +  2 


•  Angular  diameters  are  pven  in  milharcseconds  Errors  correspond  lo 
1  o 


450.0  nm  ( 1 2  mas)  and  at  533.0  nm  (1 1  mas),  while  at  656.5  nm 
it  was  essentially  unresolved  ( £  5  mas)  The  error  bars  for  April 
measurements  were  relatively  large  (±4  mas)  and  reflect  the 
relatively  low  signal-to-noise  ratio  in  the  power  spectra  due  to 
poor  seeing. 

In  May-June  we  measured  the  SN  diameter  at  four  wave¬ 
lengths:  450.0,  533.0,  656.5,  and  775  0  nm.  The  smallest  dia¬ 
meter  was  measured  at  656.5  nm  (8  +  1  mas)  The  diameters 
measured  at  450.0,  533.0,  and  775.0  nm  were  23  ±  4,  18  +  2, 
and  15+1  mas,  respectively. 

In  1987  November  diameter  measurements  were  obtained  at 
six  different  wavelengths:  530.0,  640.0,  656.6,  658.5,  700.0.  and 
850.0  nm.  The  diameters  ranged  from  15  mas  at  640.0  nm  to  20 
mas  at  533.0  nm.  The  signal-to-noise  ratios  in  these  data  were 
higher  than  in  the  data  from  the  two  previous  observing  runs. 
This  resulted  in  an  improvement  of  the  measurement  accuracy 
(making  it  ±  1-2  mas.) 

We  observed  a  substantial  asymmetry  in  the  data  recorded 
m  1988  February-March  and  1988  April  at  several  wave¬ 
lengths  (Karovska  et  al.  1988).  The  images  were  elongated 
30%-40%  along  an  axis  tilted  20c-30°  from  the  north  Images 
of  comparison  stars  and  a  Cen  A  did  not  show  this  asymmetry. 
These  results  will  be  more  extensively  described  in  a  separate 
paper. 

In  Table  2,  we  present  our  estimates  of  an  “effective”  dia¬ 
meter  for  the  SN  in  1988  February-March  and  1988  Apul 
obtained  by  fitting  the  azimuthally  averaged  power  spectra  to 
visibility  functions  computed  for  different  stellar  disks  with 
uniform  brightness  distribution 

Figures  4-6  show  examples  of  the  data  and  fitted  curves  In 
Figure  4,  we  show  the  supernova  in  Ha  from  the  1987  April 
run,  an  unresolved  comparison  star,  and  the  fitted  curves  The 
fit  is  a  horizontal  straight  line  for  both  of  them,  which  is  the 
expected  result  for  an  unresolved  star.  In  Figure  5,  we  plot  the 


Fig  4 — Power  spectra  and  best  fits  for  SN  1987A  (1987  Apnl)  and  an  unresolved  comparison  star  recorded  in  Ha  Vertical  dotted  lines  show  range  of  fitted 
region 


Fig  5— Power  spectra  and  best  fits  for  SN  1987a  Ha  data  from  1987  May-June.  198?  November.  1988  February-March.  and  1988  Apnl  observing  runs 
Vertical  dotted  lines  show  range  of  fitted  region 
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Fig  6 — Power  speura  for  SN  1987a  5J3  nm  daia  from  1988  February-March  and  Apnl  and  corresponding  btsi  fm  Vertical  lines  show  the  range  of  fitted 
region 


Ha  data  from  1987  June,  1987  November,  1988  February- 
March,  and  1988  April  and  their  fits.  The  drop  in  high- 
frequency  power  between  data  sets  corresponds  to  the 
increasing  size  in  the  SN  diameter  In  Figure  6,  we  show  the 
533  nm  data  and  fits  from  1988  February-March  and  1988 
April  The  signal-to-noise  ratio  for  these  data  was  extremely 
good  and  allowed  us  to  determine  the  SN  diameter  to  an  accu¬ 
racy  of  better  than  1  mas. 

v.  DISCUSSION 

Our  measurements  of  the  SN  diameter  on  April  2  (38  days 
after  the  explosion)  and  May-June  (95-98  days  after  the 
explosion)  show  a  substantial  dependence  on  the  wavelength  of 
observation  At  both  epochs  we  obtained  the  smallest  diameter 
at  656.5  nm  and  the  largest  at  45C.O  and  533.0  nm.  The  656.5 
nm  filter  was  located  redward  from  the  Ha  absorption 
minimum  and  encompassed  a  large  fraction  of  the  P  Cygni  Ha 
emission  feature  (Phillips  1987).  The  450.0  nm  filter  is  in  the 
absorption  line  corresponding  to  Ba  II  (Phillips  1987),  while  no 
prominent  absorption  or  emission  line  could  be  distinguished 
in  the  533  nm  region.  In  May-June  we  observed  the  SN  using 
an  additional  filter  centered  on  the  O  i  absorption  feature  at 
775  0  nm.  At  this  wavelength,  the  diameter  is  smaller  than  the 
one  measured  at  450.0  and  533  0  nm,  but  almost  twice  as  large 
as  the  diameter  measured  at  656.5  nm. 

The  wavelength  dependence  of  the  measured  diameters  is 
probably  due  to  the  stratification  of  the  expanding  envelope  of 
the  SN.  Evidence  of  stratification  in  the  SN  shell  has  already 
been  found  from  spectral  observations  of  the  absorption 
minima  of  several  spectral  lines  (Ha,  H/?,  Hy.  Ca  u,  Na  i)  during 
the  first  month  of  expansion  (Hanuschik  and  Dachs  1987a,  b ). 
However,  the  (act  that  the  angular  diameter  measured  in  the 
O  i  absorption  trough  is  twice  that  measured  in  the  Ha  emis¬ 


sion  line  requires  that  hydrogen  lie  within  the  oxygen  shell,  a 
result  that  is  not  predicted  by  current  SN  atmosphere  models. 

It  is  interesting  to  note,  however,  that  the  speckle  diameters 
measured  at  the  656.5  nm  spectral  bandpass  in  1987  April  and 
May-June  are  consistent  with  the  estimates  of  the  effective 
diameters  of  the  line-forming  regions  at  both  epochs  Figure  7 


Fic  7 — Diameters  vs  lime  for  SN  1 987 A  calculated  using  velocities  esti¬ 
mated  from  minima  of  a  various  spectral  lines  and  the  blackbody  photometric 
diameter  estimates  for  the  first  130  davs  since  explosion  (Phillips  ei  al  1988) 
The  two  measurements  in  Ha  and  at  533  nm  are  also  shown 
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shows  the  estimates  of  the  effective  diameters  of  the  line- 
forming  regions  which  were  obtained  by  integrating  the  radial 
velocities  corresponding  to  the  absorption  minima  of  several 
spectral  lines  during  the  first  130  days  after  the  explosion 
(Hanuschik  and  Dachs  1987a,  b ;  Phillips  ei  al.  1988).  On  the 
38th  day  after  the  explosion,  the  estimated  angular  diameter  of 
the  Ha  she'll  (i.e.,  the  diameter  of  the  region  where  Ha  becomes 
optically  thick)  is  6  mas.  and  we  measured  an  angular  diameter 
<,5  mas  at  656.5  nm.  The  angular  diameter  corresponding  to 
the  H/I  and  Hy  “shell"  is  around  5  mas.  The  estimated  dia¬ 
meters  derived  from  the  radial  velocity  of  the  minimum  of  the 
Ba  r  absorption  line  is  1  R  mas  On  the  98th  day  after  the 
explosion  we  obtained  an  1  nas  diameter  at  656.5  nm,  which  is 
in  close  agreement  with  the  estimated  effective  diameter  of  9 
mas  for  the  H/i  and  Hy  “shells,”  and  a  7  mas  diameter  for  the 
Ba  u 

The  diameters  measured  in  the  continuum  bandpasses  in 
1987  April  and  1987  May-June  are  substantially  larger  than 
those  obtained  from  photometry  and  a  simple  blackbody  fit  to 
the  observed  spectrum  (Catchpole  el  al.  1987;  Danzinger  ei  al. 
1987,  Whitelock  ei  al.  1988)  These  estimates  of  the  SN  dia¬ 
meter  are  shown  in  Figure  7.  Larger  continuum  diameters 
obtained  using  speckle  interferometry  may  be  a  result  of  elec¬ 
tron  scattering  in  the  SN  envelope  Lucy  (1987)  computed  the 
"last-scattering"  diameters  (scaled  to  the  blackbody  photo- 
spheric  diameter)  for  the  38th  and  98th  day  after  explosion. 
’  hese  computations  were  performed  for  the  appropriate  band- 
passes  of  the  speckle  observations.  Calculated  diameters  are 
2-4  times  largei  than  the  blackbody  diameter  estimates,  but 
they  are  still  1 .5—3  times  smaller  than  the  diameters  obtained 
using  speckle  interferometry  This  suggests  that  the  phenome¬ 
non  of  electron  scattering  cannot  entirely  account  for  the 
larger  size  of  the  SN  measured  by  speckle  interferometry. 
These  models  assume  a  uniform  distribution  of  material  in 
spherically  symmetric  shells.  Clumpiness  or  asymmetries  in  the 
distribution  of  the  material  may  well  produce  very  different 
results 


The  observations  in  1987  November,  1988  February- 
March,  and  1988  April  were  obtained  during  the  period  when 
the  SN  was  already  in  its  "nebular”  phase.  The  SN  diameter 
measured  at  656.5  nm  showed  clear  temporal  evolution  when 
compared  with  our  earlier  measurements  at  the  same  wave¬ 
length.  Figure  8  shows  the  evolution  of  the  Ha  diameter  as 
measured  by  speckle  interferometry  during  the  first  41!  days 
after  the  explosion.  In  addition  to  our  measurements  we  show 
the  speckle  measurements  near  the  center  of  the  Ha  line 
obtained  by  Wood  er  al.  (1988)  on  1987  December  12  and  13, 
and  1988  February  20  using  the  Anglo-Australian  Observatory 
4  m  telescope.  The  uniform  disk  angular  diameters  estimated 
for  these  two  epochs  are,  respectively.  23.1  ±  1.6  and 
23.9  +  0.9  mas.  We  performed  a  linear  least-squares  fit  to  these 
data  assuming  that  the  diameter  at  the  moment  of  the  explo¬ 
sion  was  zero  If  the  increase  of  the  diameter  reflects  the  expan¬ 
sion  of  the  SN  atmosphere,  then  the  slope  of  the  fitted  straight 
line  gives  a  mean  velocity  of  expansion  of  approximately  2850 
km  s'1.  Our  measurement  from  February-March  seems  to 
underestimate  the  diameter  of  the  SN  at  that  epoch.  This  mea¬ 
surement  has  been  affected  by  the  asymmetry  in  those  data 

Diameter  measurements  from  the  last  three  observing  runs 
still  show  some  dependence  on  the  wavelength  of  observation, 
though  the  differences  between  the  diameters  measured  in  Ha 
and  in  the  other  wavelengths  are  not  as  dramatic  as  they  were 
in  1987  April  and  May-June  Diameter  measurements  in  the 
S33  nm  continuum  from  1987  May-June  and  November  and 
1988  February  are  practically  indistinguishable  from  one 
another.  This  is  also  true  for  the  diameters  measured  at  640, 
700,  and  850  nm  in  1987  November  and  1988  February- 
March  If  the  533  nm  measuiements  correspond  io  the  size  of 
the  SN  photosphere  then  we  conclude  that  they  do  not  indicate 
any  recession  at  these  epochs 

At  this  point,  it  is  difficult  to  find  a  simple  interpretation  for 
all  our  results.  New,  more  detailed  models  for  the  SN  arc 
required  We  plan  to  continue  to  monitor  the  changes  in  the 
SN  diameter  over  a  large  spectral  range  and  to  follow  the 
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Fic  8  —  E  volution  of  the  SN  diameter  during  the  first  41 1  days  after  the  explosion  as  observed  by  speckle  interferometry  techniques  wilh  J  o  errurs  and  best-fit 
straight  line  giving  a  velocity  of  2850  km  s'1 
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evolution  of  the  expanding  shell.  If  the  present  slow  decline  in 
brightness  of  the  SN  continues,  we  can  observe  the  SN  for 
several  years  belore  the  magnitude  limit  of  the  speckle  process 
is  reached. 
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ABSTRACT.  The  speckle  imaging  program  at  the  CfA  has  developed  a  set  of  tools  for  the 
implementation  of  speckle  image  reconstruction  and  has  applied  these  tools  to  a  wide  range  of  scientific 
programs.  Two  key  elements  in  this  program  have  been  the  development  of  a  2-dimensional  photon 
counting  camera,  the  PAPA  detector,  and  the  implementation  of  image  reconstruction  procedures  using 
the  Knox-Thompson  algorithm.  This  paper  discusses  the  current  status  of  the  PAPA  detector  and  the 
details  of  our  data  processing  methods  and  numerical  algorithms. 


1.  Introduction 

Speckle  interferometry  is  proving  to  be  an  invaluable  tool  for  astronomical 
research,  allowing  measurements  of  a  wide  range  of  scientifically  interesting  objects 
that  can  be  made  in  no  other  way.  Most  of  die  scientific  results  from  speckle  to  date  are 
in  the  areas  of  binary  star  orbits  and  stellar  angular  diameters.  These  objects  are 
generally  bright  and  have  relatively  simple  geometries.  Extending  these  techniques  to  a 
wider  range  of  problems  entails  the  observation  of  much  fainter  sources  with,  in  many 
cases,  more  complex  geometries.  At  the  Center  for  Astrophysics  (CfA),  we  have  been 
concentrating  on  the  development  of  detectors  and  algorithms  that  allow  application  of 
the  speckle  process  to  image  reconstruction  of  faint,  extended  sources.  The  Knox- 
Thompson  (K-T)  algorithm  (Knox,  1976)  is  one  of  a  class  of  approaches  to  preserving 
the  phase  in  the  object  Fourier  transform,  allowing  true  image  reconstruction  as  an 
output  of  the  process.  The  CfA  application  of  the  K-T  algorithm  has  proven  to  be  far 
more  robust  and  accurate  in  its  phase  recovery  process  than  Knox's  original 
implementation  due  to  several  improvements  in  the  technique.  These  improvements 
include  the  averaging  of  complex  vector  differences  rather  than  phase  differences  and 
application  of  a  least  squares  error  distribution  in  the  phase  estimation  process. 

A  crucial  technology  that  allows  the  application  of  K-T  to  faint  astronomical 
sources  has  been  the  development  of  the  PAPA  2-dimensional  photon  counting  camera. 
This  detector  records  the  spatial  coordinates  and  time-of-arrivad  of  each  detected  event, 
allowing  the  extremely  accurate  corrections  for  photon  noise  biases  required  for 
accurate  reconstruction  from  the  K-T  process.  The  PAPA  has  proved  to  be  a  powerful 
and  essential  tool  in  the  application  of  speckle  imaging,  providing  the  crucial  linearity 
and  dynamic  range  needed  for  accurate  measurements  of  source  characteristics.  In  this 
paper,  the  details  of  the  CfA  implementation  of  speckle  imaging  using  the  PAPA 
derector  and  modified  versions  of  the  K-T  algorithm  are  discussed. 


2.  Data  Recording 
2.1  TOE  PAPA  DETECTOR 


,oot.  PAPA  (Precision  Analog  Photon  Address)  detector  (Papaliolios,  et  al, 
1985)  allows  recording  of  the  address  (position)  and  time  of  arrival  of  each  detected 
photon.  A  schematic  diagram  of  the  PAPA  is  shown  in  figure  1.  The  front  end  of  the 
camera  is  a  high  gain  image  intensifier  which  produces  a  bright  spot  on  its  output 
phosphor  for  events  detected  by  the  photocathode.  The  back  face  (phosphor)  of  the 
intensifier  is  then  reimaged  by  an  optical  system  which  is  made  up  of  a  large  collimating 
lens  and  an  array  of  smaller  lenses.  Each  of  the  small  lenses  produces  a  separate  image 
of  the  phosphor  on  a  binary  mask.  Behind  each  mask  is  a  field  lens  which  relays  the 
pupu  of  the  small  lens  onto  a  small  photomultiplier  (PMT) 
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Figure  1.  The  PAPA  Detector  and  Gray  Scale  Masks 


The  plan  for  the  masks  is  shown  at  the  bottom  of  figure  1.  Each  mask  gives 
one  bit  of  the  digital  address  of  the  detected  photon  either  in  the  x  or  y  directions.  One 
mask  is  completely  clear,  and  this  channel  is  called  the  strobe  channel.  When  the  strobe 
sees  a  pulse,  signaling  the  detection  of  an  event,  the  other  channels  are  switched  on  and 
interrogated.  A  discriminator  level  is  set  for  all  the  other  channels,  scaled  to  the  height 
of  the  pulse  detected  by  the  strobe  channel.  If  the  photon  image  falls  on  the  opaque  part 
of  the  mask  in  a  given  channel,  the  output  of  the  PMT  will  be  below  the  discriminator 
level  and  that  bit  is  set  to  zero.  If  the  photon  image  falls  on  a  clear  area,  the  PMT  sees  a 
pulse  above  the  threshold  and  the  address  bit  is  unity  .  The  masks  use  a  Gray  (instead 
of  binary)  code  which  ensures  that  the  mask  stripes  do  not  have  edges  located  in  the 
same  place  in  the  field. 

The  CfA  PAPA  has  19  channels  giving  9  bit  by  9  bit  resolution  (512x512)  plus 
the  strobe  channel.  We  can  currently  record  up  to  100,000  event  addresses  per  second 
from  the  camera.  The  inherent  limit  on  the  PAPA  data  rate  is  set  by  the  decay  time  of 
the  intensifier  phosphor,  (about  200  ns  for  P47),  which  suggests  that  the  camera  could 
operate  at  a  rate  over  1  million  counts  per  second.  The  current  camera  uses  only  a 
single  Gen  II  (channel  plate)  intensifier,  which  has  only  3  to  4%  quantum  efficiency. 
New  tubes  which  couple  a  Gen  I  diode  in  front  of  the  Gen  II  have  recently  improved 
the  peak  quantum  efficiency  to  over  10%.  The  high  data  rate  requires  special  recording 
capabilities.  At  CfA,  we  have  built  a  VCR  digital  recording  system  which  encodes  the 
digital  output  of  the  camera  onto  a  video  carrier,  allowing  the  storage  of  two  hours  of 
telescope  data  on  a  single  VHS  tape  (Ebstein,  1986).  The  data  must  be  stripped  off  the 
tape  and  reconstituted  into  a  conventional  digital  form  for  numerical  processing. 

The  PAPA  is  a  highly  linear  detector  with  a  very  large  dynamic  range.  Unlike 
other  2-D  photon  detectors,  the  data  rate  and  resolution  are  independent  of  the  image 
intensity  distribution.  The  major  problem  with  the  PAPA  is  that  aberrations  in  the 
optics  cause  errors  in  the  detected  photon  positions  which  are  dependent  on  the  position 
in  the  field.  This  shows  up  as  a  weak  tartan-like  pattern  in  the  background.  Recording 
of  data  with  uniform  illumination  allows  flat  field  corrections  of  the  data  that  essentially 
eliminate  all  traces  of  the  baci'  "round  errors.  Details  of  the  flat  fielding  procedure  are 
given  in  section  2.3. 


2.2  RECORDING  REQUIREMENTS 

Stellar  speckle  data  recorded  to  be  reduced  using  K-T  have,  in  almost  all 
respects,  the  same  requirements  as  those  for  standard  speckle  interferometry.  The  key 
requirements  are 

-  Short  exposure  times  which  "freeze"  the  atmospheric  turbulence.  With  the  PAPA 
detector,  the  time  of  arrival  of  each  event  is  recorded,  so  photons  may  be  grouped  into 
frames  in  a  way  which  maximizes  the  signal-to-noise  ratio  in  the  integrated  power 
spectrum.  Tests  are  usually  run  on  comparison  star  data  at  a  series  of  "exposure” 
times,  and  the  signal-to-noise  calculated  at  several  frequencies  in  the  power  spectrum. 

-  Diffraction  Sampling.  Magnifying  optics  are  used  in  front  of  the  camera  to  rescale  the 
image  so  that  there  are  at  least  2  (and  preferably  3)  pixels  for  each  telescope  diffraction 
element. 


-  Dispersion  Correction.  Counter-rotating  computer  controlled  dispersion  correcting 
prisms  are  used  to  compensate  for  the  atmospherically  induced  dispersion  at  zenith 
angles  larger  than  a  few  degrees. 

-  Isoplanatism.  The  processed  field  size  must  be  small  enough  so  that  the  atmosphere- 
telescope  point  spread  function  is  invariant  over  the  entire  field. 

-  Spectral  Bandpass.  The  bandpass  must  be  narrow  enough  to  provide  temporal 
coherence  over  the  whole  image  plane.  A  good  rule  of  thumb  requirement  is  that 

Bandpass  =  1000A/(S-  D)  where  S  is  the  seeing  in  arcseconds  and  D  is  the  telescope 
diameter  in  meters. 

-  Seeing.  One  requirement  for  Knox-Thompson  th'U  is  not  required  for  speckle 
interferometry  is  that  the  diameter  of  the  seeing  disk  be  less  than  1/2  the  field  size,  since 
this  sets  the  scale  of  phase  differences  in  the  K-T  implementation.  In  addition,  the 
number  of  frames  required  for  a  given  signal-to-noise  ratio  increases  as  the  4th  power 
of  the  seeing,  so  good  seeing  is  imperative  for  good  results. 

-  Photon  noise  limited  detection.  Correction  for  the  photon  noise  bias  is  a  crucial  step 
in  the  implementation  of  K-T.  The  PAPA  detector  and  other  detectors  like  it  are  critical 
in  allowing  exact  estimation  and  correction  of  these  biases,  unlike  other  2-D  analog 
cameras. 


2.3  FLAT  FIELDING 

The  procedure  for  flat  fielding  the  data  is  to  first  generate  a  flat  field  "mask" 
using  a  combination  of  "white-field"  data  (in  which  the  camera  is  looking  at  the  dome  or 
an  internal  uniform  source)  and  a  filtered  version  of  the  long-exposure  image  of  the 
object  being  flat  fielded  (Ebstein,  1987a).  Use  of  the  object  data  itself  allows  for  any 
dynamic  changes  in  the  flat  field  and  reduces  the  requirement  of  recording  "white-field" 
data  to  about  once  per  night.  The  synthesized  fiat  field  is  then  applied  to  each  frame, 
after  the  photon  addresses  have  been  binned  into  individual  frames.  One  ends  up  with 
frames  in  which  each  photon  is  given  a  fractional  height  proportional  to  the  flat  field 
amplitude  at  that  position.  The  steps  in  the  operation  are 

1 .  Integrate  the  "white  field"  data  into  a  long  exposure  and  low  pass  filter  with  a  cutoff 
frequency  of  X  /  Tq  to  eliminate  the  high  frequencies  in  this  image, 

2.  Generate  the  long  exposure  of  the  object  using  direct  integration  of  the  photon  data. 

3.  Divide  the  object  long  exposure  (step  2)  by  a  low  pass  filtered  (cutoff  =  X  /  Tq  ) 

version  of  the  object  long  exposure.  This  operation  eliminates  the  object  related  low 
frequencies  from  the  result. 

4.  Multiply  the  white  field  image  by  the  result  of  the  division  performed  in  step  3. 

5.  Take  the  inverse  of  the  result  of  step  4  and  then  apply  a  border  mask  which  rolls  off 
the  region  near  the  edges  where  there  is  a  low  signal-to-noise  ratio. 


6.  In  composing  frames  from  the  photon  address  data,  replace  each  photon  with  the 
value  of  the  flat  field  at  that  address. 

The  corrected  frames  then  have  the  form 


I(r)  =  2  FF(rj)  8(  r-rj ) 
.  where  FF(fj)  is  the  flat  field. 


3.  The  Knox-Thompson  Algorithm 
3.1  SPECKLE  INTERFEROMETRY 

Short  exposure  images  recorded  in  the  image  plane  of  a  large  telescope  have  the  form 

I(r)  =  0(r)®  S(r)  3.1 

where  I(r)  is  the  image,  O(r)  is  the  object  intensity  distribution  and  S(r)  is  the 
telescope-atmosphere  point  spread  function.  ®  represents  the  convolution  operation. 

Fourier  transforming 

1(f)  =  6(f)  •  S(f)  3.2 


Here  O(f)  is  the  object  spectrum  and  S(f)  is  the  transfer  function.  In  conventional 
speckle  interferometry,  the  ensemble  averaged  power  spectrum  is  obtained  for  a  large 
set  of  short  exposure  images 

<II(f)|2>=  10(f)!2  •  <IS(f)l2>  3.3 

Here  <  >  indicates  the  ensemble  average  and  1 1  the  modulus. 


3.2  KNOX-THOMPSON  IN  FOURIER  SPACE 
In  the  Knox-Thompson  process  we  save  three  quantities 

<A(f)>  =  <  11(f)!2  >  -  <2  FF2(fj)>  3.4 

j  J 

<Du(f)>  =  <I*(f)  I(f+Afu)  > 

=  <  II(f)l2>  •  <  exp{-i[<J>(f)  -  0(f  +  Afu)]J  > 

-  <2  FF2(fj)  •  exp(i[Afu  fj]> 


3.5 


3.6 


<Dv(f)>  =  <  I*(f)  I(f+Afv)  > 

=  <  II(f)l2>  <  cxp(-i[$(D  -  W  +  Afv)])  > 

-  <Z  FF2(f:)  •  exp(i[Afv  f:]> 

J  J  J 

The  quantity  in  equation  3.4  is  die  conventional  result  of  speckle  interferometry, 
photon  noise  bias  corrected  power  spectrum.  After  ensemble  averaging,  the  power 
spectrum  may  be  Fourier  transformed  back  to  image  space,  producing  an 
autocorrelation  with  centro-symmetry.  While  information  such  as  binary  separations, 
position  angles  and  relative  intensities  of  the  components  can  be  extracted  from  the 
autocorrelation,  in  many  cases  important  information  is  lost,  particularly  for  complex 
extended  objects.  The  two  additional  terms  saved  in  the  K-T  process  are  given  in 

equations  3.5  and  3.6.  They  represent  the  complex  autocorrelation  where  Afu  and  Afv 
are  small  (compared  to  the  atmospheric  correlation  length,  tq)  constant  displacements  in 
Fourier  space. 


Figure  2.  Speckle  Transfer  Function  +  Photon  Noise  Bias 


Two  crucial  steps  in  the  implementation  of  K-T  are  the  precise  correction  of  the 
photon  noise  biases  and  the  correction  for  the  speckle  transfer  function.  The  reasons 
for  this  are  demonstrated  in  figure  2.  For  faint  objects,  the  relatively  small  signal  sits 
on  a  huge  bias  which  must  be  subtracted  for  an  accurate  reconstruction.  If  the  photons 
are  not  delta  functions,  the  bias  will  have  a  shape  which  must  be  accurately  known. 
Equivalent  bias  adds  large  phase  errors  to  the  cross  spectrum  terms  and  severely 
degrades  the  reconstructed  phases.  Since  most  of  the  power  in  the  speckle  power 

spectrum  is  located  at  frequencies  below  X/rg  ,  the  long  exposure  cutoff,  accurate 

estimation  of  the  transfer  function  must  be  made  from  comparison  star  data.  Small 
changes  in  rg  between  the  recording  of  the  object  and  the  comparison  star  will  yield 

large  errors  in  amplitude  calibration. 

The  photon  noise  bias  is  easily  corrected  for  data  from  a  photon  counting 
camera  such  as  the  PAPA,  since  the  photon  count  is  exactly  known,  and  the  centroid  of 
each  photon  image  is  detected  as  a  delta  function  (unlike  analog  cameras  where  a 
photon  image  spreads  over  several  pixels  and  the  photon  count  may  only  be  estimated). 
When  the  frames  are  flat  fielded  and  the  photons  are  no  longer  unit  amplitude,  the  bias 
correction  terms  are  derived  by  first  integrating  the  square  of  the  long  exposure  image. 

The  Fourier  transform  of  this  image  is  then  evaluated  at  f  =  0,  f  »  Afu  and  f  ■  Afv  . 
These  quantities  are  then  subtracted  from  the  con-esponding  integrated  power  and  cross 

A  A 

spectra,  as  indicated  in  equations  3.3  to  3.5.  Du  and  Dv  may  be  expanded.  For 
example 


<Du(f)>  *  0(f)  •  0(f+Afu)exp{i[  0(f)  -  <|>(f-Afu)]} 

x  <  T(f)  •  T(f+Afu)  >  •  <  exp{i[  ^(f)  -  T(f-Afu)j)  >  3.7 


A  A  A  > 

But  <T(f)  •  T(f+Afu)  >  -  <  IT(f)l^>  which  is  just  the  speckle  interferometry  transfer 
function,  a  quantity  which  may  be  estimated  from  an  unresolved  star  close  in  angular 
position  to  the  object,  and  PF(f)  -  'F(f-Afu)]  is  the  atmospheric  phase  difference, 

which  is  a  mean  zero  stochastic  process,  if  the  telescope  adds  no  asymmetric 
aberrations.  Therefore,  after  averaging  we  are  left  with  two  complex  arrays  which  are 
weighted  by  the  speckle  transfer  function  and  which  contain  phase  difference 
information  in  the  u  and  y  directions.  We  will  use  the  power  spectrum  data  to  extract 
the  amplitudes  of  the  image  Fourier  transform  and  these  phase  difference  arrays  to 
recover  the  image  phase. 

A  A 

After  averaging,  we  extract  phases  from  the  Du  and  Dv  arrays  using  an  iterative 

least  squares  approach.  The  details  of  the  phase  recovery  process  are  discussed  in 
section  3.5 


3.3  KNOX-THOMPSON  IN  IMAGE  SPACE 

Data  recorded  in  a  photon  counting  mode  may  easily  be  integrated  in  image 
space  using  the  photon  coordinates.  This  is  advantageous  for  the  implementation  of 
real-time  correlators  and  even  for  some  computation  schemes.  The  computation  time 


required  to  perform  the  image  plane  correlation  increases  as  the  square  of  the  number  of 
photons,  while  the  Fourier  transform  is  performed  frame  by  frame,  even  when  the 
frame  is  nearly  empty.  Therefore,  depending  on  the  processor,  the  image  plane 
operation  will  be  more  efficient  for  small  numbers  of  photons/frame  (usually  on  the 
order  of  a  few  hundred  to  a  thousand)  while  an  FFT  will  be  better  at  higher  light  levels. 

Again,  three  quantities  are  ensemble  averaged,  this  time  in  image  space 


<A(r)>  =  <  I  6(rn  -rn-  -  r)  > 

3.8 

nn 

<K  (r)>  =  <  I,  S(rn  -  rn*  -  r)  •  exp(-iAurn0  > 

3.9 

nn 

<Kv(r)>  =  <  I,  8(rn  -  rn>  -  r)  exp(-iAvrn0  > 

3.10 

In  this  case,  no  bias  terms  are  accumulated  if  we  simply  discard  all  the  self-correlation 
terms  (n=n’),  so  no  bias  correction  is  needed.  The  complex  arrays  may  be  saved  as  real 
and  imaginary  arrays,  requiring  only  a  look-up  table  having  n  values  (where  n  is  the 
array  dimension)  for  the  complex  coefficients. 

Fourier  transforming  these  quantities  after  the  ensemble  average  converts  them 
to  be  equivalent  to  the  Fourier  integration  results,  allowing  the  phase  and  amplitude 
reconstruction  procedures  described  in  sections  3.4  and  3.5  to  be  applied. 


3.4  AMPLITUDE  RECONSTRUCTION 

If  speckle  data  is  recorded  for  an  unresolved  star  as  well  as  for  the  object,  the 
this  comparison  star  data  may  be  used  to  correct  the  speckle  transfer  function.  In 
general,  the  comparison  star  is  chosen  to  be  as  close  in  angular  position  to  the  object  as 
possible.  Data  sets  on  the  comparison  star  are  taken  interspersed  with  data  on  the  object 
so  that  seeing  changes  are  minimal.  We  have  found  that  taking  5  minutes  of  data  on  the 
comparison  star,  then  10  minutes  on  the  object,  and  then  5  minutes  on  the  comparison 
star  again  is  about  optimal  in  terms  of  efficient  use  of  telescope  time  and  minimizing 
seeing  changes. 

In  the  data  processing,  the  comparison  star  data  is  integrated  using  the  same 
procedures  for  the  object.  The  last  operation  in  the  reconstruction  procedure  is  to  divide 
the  amplitudes  from  the  object  spectrum  by  the  amplitudes  from  the  comparison  star.  A 
similar  operation  is  used  to  correct  the  object  phases  using  the  comparison  star 
reconstructed  phases.  The  corrected  amplitudes  are  then  rolled  off  using  a  low  pass 
filter  (usually  in  the  form  of  an  "optimum"  filter  which  minimizes  sidclobes  in  the 
image). 


3.5  PHASE  RECONSTRUCTION 


The  scheme  for  phase  reconstruction  from  the  complex  cross  spectrum  is  shown 
in  figure  3.  The  procedure  for  phase  recovery  consists  of  setting  one  point  (usually  the 
f  =  0  point)  equal  to  unit  amplitude  and  zero  phase.  First  estimates  for  the  phasors 
along  the  fust  row  and  first  column  are  then  calculated  as  the  product  of  the  phasor  at 
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Figure  3.  Phase  Integration  Procedure 

the  previous  point  multiplied  by  the  appropriate  Du  or  Dv  term,  normalized  back  to  a 

phasor.  Use  of  the  complex  amplitude  in  the  phase  integration  produces  a  natural 
weighting  of  the  phase  reconstructor,  proportional  to  the  signal  along  that  path.  The 
interior  points  are  then  calculated  using  a  two  term  average.  Set  for  example  the 
calculation  for  P22  in  Figure  3.  These  operations  provide  a  fairly  crude  fust  estimate  of 


the  phase,  leaving  large  cumulative  errors  at  higher  frequencies  (further  away  from 
f=0).  This  error  is  then  redistributed  over  the  entire  frequency  plane  using  an  iterative 
algorithm  which  converges  to  a  result  equivalent  to  the  least  squares  solution  for  the 
phasors.  The  iteration  is  demonstrated  in  figure  3  to  obtain  a  new  estimate  of  P23.  The 
iteration  starts  at  the  f  =  0  point  and  works  its  way  out  to  the  maximum  frequency 
following  a  spiral  path.  Successive  clockwise  and  counter-clockwise  spiral  paths  are 
followed.  The  iteration  continues  until  a  stable  solution  is  reached,  usually  between  20 
to  100  cycles.  The  resulting  phasors  are  corrected  for  the  effects  of  asymmetric 
telescope  aberrations  by  a  division  with  the  reference  star  phasors.  The  final 
reconstructed  phasors  are  then  multiplied  by  the  deconvolved  amplitudes  and  Fourier 
transformed  back  to  image  space  to  produce  an  image. 


3.6  EXTENDED  KNOX  THOMPSON 

Improved  signal-to-noise  can  be  obtained  with  K-T  by  calculating  additional 
complex  cross  spectrum  arrays  for  other  scales  of  Afx  and  Afv.  These  phasor 

estimates  may  then  be  averaged  with  the  original  estimates,  using  a  weighting 
proportional  to  the  degree  of  atmospheric  correlation  at  those  spacings.  In  general,  only 
a  few  additional  terms  contribute  to  improved  signal-to-noise,  since  except  under 
extremely  good  conditions,  the  atmospheric  correlation  drops  rapidly  to  zero  at  larger 
separations.  The  diagonal  terms  add  most  of  the  additional  information. 


3.7  CLEANUP  ALGORITHMS 

Several  algorithms  are  available  for  "cleaning  up"  the  images  after  the  K-T 
reconstruction.  These  include  Fienup,  Constrained  Iterative  Algorithm  (CIA),  Iterative 
Deconvolution,  and  Maximum  Entropy.  Some  success  has  been  obtained  with  all  of 
them,  however  they  are  all  non-linear,  so  while  they  may  improve  the  image 
morphology,  they  do  not  appear  to  preserve  the  accuracy  of  the  amplitudes.  The 
original  reconstructions  (and  for  simple  objects,  the  autocorrelation  images)  seem  to 
give  the  most  accurate  amplitude  estimates.  A  summary  of  the  algorithms  and  their  uses 
follows: 

-  The  Fienup  algorithm  (Fienup,  1978)  applies  image  plane  constraints  including 
positivity  and  finite  support  to  adjust  the  phases  while  keeping  the  original  amplitudes. 
This  algorithm  was  originally  suggested  for  obtaining  phases  when  only  the  amplitudes 
are  known.  Starting  with  the  K-T  phase  estimates  and  high  signal-to-noise 
deconvolved  amplitudes,  this  technique  usually  converges  in  a  few  cycles,  reducing 
background  with  little  effect  on  the  photometric  accuracy.  The  success  of  the  procedure 
is  highly  dependant  on  the  signal-to-noise  in  the  amplitude  spectrum. 

-  Constrained  Iterative  Algorithm  (CIA)  is  designed  to  clean  up  the  amplitudes  by 
applying  a  set  of  image  and  Fourier  plane  constraints  (Ebstein,  1987b).  Iteration 
between  Fourier  space  and  Autocorrelation  space  constrains  the  deconvolution  to  be 
between  signal-to-noise  bounds,  and  the  autocorrelation  to  be  positive  and  to  have  a 
finite  support.  This  algorithm  is  effective  for  improving  noisy  amplitudes,  but  it  can 
produce  spurious  peaks  in  the  image  and  docs  not  preserve  amplitude  ratios  with  great 
accuracy. 


-  Iterative  Deconvolution  (Ayers  and  Dainty  1988)  iterates  between  new  estimates  of 
the  image  and  of  the  point  spread  function  using  similar  image  plane  constraints  as 
Fienup.  This  approach  appears  to  produce  very  clean  deconvolutions,  with  no  side 
band  problems.  Its  linearity  is  untested  at  this  time. 

-  Maximum  Entropy  (MEM)  Many  different  implementations  of  maximum  entropy 
have  been  introduced.  This  technique  has  proven  to  be  very  successful  for  improving 
radio  maps.  Again  this  technique  has  been  reasonably  successful  at  cleaning  up  when 
good  estimates  of  the  images  are  provided  from  K-T.  However,  accurate  amplitude 
ratios  do  not  appear  to  be  preserved. 


3.8  CONVERGENCE 

The  convergence  of  the  phase  from  the  K-T  integration  may  be  calculated 
(Nisenson  and  Papaliolios,  1983)  in  a  manner  very  similar  to  the  calculation  for  the 
amplitudes  (Goodman  and  Belsher,  1976).  The  approach  is  to  calculate  the  variance  of 

A 

the  complex  autocorrelation.  For  example,  the  variance  for  Du(f)  is  given  by 


o2(0  =  <  IDu(f)  |2  >  -  kDu(f)  >|2 


3.11 


From  this  analysis,  we  can  show  that  the  number  of  frames  required  for 
convergence,  M,  is  related  to  the  number  of  detected  photons,  P,  the  number  of 

A 

speckles  in  the  seeing  disk,  ns,  the  object  spectrum,  O  ,  and  the  allowed  phase  error,  E, 
at  the  spatial  frequency  f/fmax  where  fmax  is  the  frequency  cutoff  (usually  the 
telescope  diffraction  limit).  Then 
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A  reasonable  requirement  for  convergence  is  to  choose  e  to  be  1/10  X  at  a  frequency  of 
0.9  (90%  to  the  diffraction  limit),  then 

M=125[~-r2  3.13 

where  ns  -  2  •  (D/rg)  2 

We  point  out  that  the  number  of  frames  required  for  convergence  increases  as 
the  4th  power  of  the  seeing  (rg).  It  is  obvious  that  the  magnitude  limit  of  the  process  is 

highly  dependent  on  the  quality  of  the  seeing. 

As  an  example,  with  a  13th  magnitude  star  we  detect  10000  photons  per  second 
with  a  4-meter  telescope,  an  overall  detection  efficiency  of  5%  (optical  and  quantum), 
and  1  arcsecond  seeing  (rg  =  10  cm).  Then  we  will  have  about  3700  speckles  and  we 

find  for  M  from  equation  3.15:  M  =  170,000  frames  or  1/2  hour  observation  with  a  full 
duty  cycle  camera.  Of  course,  even  for  faint  objects,  substantial  enhancement  can  be 


obtained  at  lower  spatial  frequencies.  The  real  criteria  are  set  by  the  specific  scientific 
goals  of  a  given  program. 


3.9  KNOX-THOMPSON  MISCONCEPTIONS 

The  original  implementation  of  K-T  (Knox,  1976)  had  many  problems  when  it 
was  applied  to  real  noisy  data.  Kno>  averaged  phases  differences  (rather  than  vector 

differences).  This  many  times  re:  dte  i  in  2n  ambiguities  in  the  average.  Zeroes  in  the 
object  transform  yield  indeterminate  phases  and  phase  differences  and  a  breakdown  in 
the  phase  integration  process.  In  addition,  the  step  of  integrating  from  phase 
differences  to  phases  was  performed  by  a  simple  summing  process  from  low  to  high 
frequencies,  resulting  in  an  integrated  error  at  the  higher  frequencies.  Finally,  Knox 
believed  that  centroiding  was  necessary  for  each  image  before  integration,  adding  to  the 
computation  time  and  increasing  the  error  for  faint  objects  due  to  poor  centroiding 
caused  by  photon  noise. 

In  the  CfA  implementation  of  Knox-Thompson,  the  problems  of  2n  ambiguities 
and  transform  zeroes  are  eliminated  by  averaging  and  integrating  complex  vector 
differences,  rather  than  phase  differences.  The  iterative  least  squares  integration 
described  in  section  3.5  distributes  the  errors  due  to  noise  over  the  entire  frequency 
plane.  It  averages  over  all  phase  paths  with  a  weighting  proportional  to  the  integrated 
amplitude  along  that  path.  In  general,  the  reconstructed  phase  will  have  an  even  better 
accuracy  than  the  vector  differences,  since  many  different  estimates  along  many 
independent  (or  partially  independent)  paths  are  averaged,  reducing  the  error.  It  should 
be  pointed  out  that  this  iterative  technique  can  also  be  used  for  one-dimensional  data, 
though  it  will  only  distribute  the  enror  over  all  frequencies  (in  the  least  squares  sense), 
not  reduce  the  total  error,  since  in  this  case  there  is  only  one  path  for  the  integration. 
Finally,  it  is  easy  to  show  that  poor  centroiding  results  in  a  scalar  error  in  the  vector 
differences  which  averages  independently.  This  error  may  be  accurately  estimated  from 
the  direct  sum  of  the  input  frames  and  then  removed  from  the  vector  differences  in  the 
ensemble,  before  the  phase  integration  step. 

One  other  misconception  about  K-T  is  that  it  is  sensitive  to  telescope 
aberrations.  It  has  been  shown  (Barakat  and  Nisenson,  1981)  that  K-T  is  only 
sensitive  to  odd  order  aberrations  (e.g.  coma)  and  not  to  defocussing,  astigmatism, 
spherical,  etc.,  and  that  the  residual  effects  of  the  aberrations  may  be  calibrated  out 
using  comparison  star  data  recorded  through  the  same  optical  system. 


4.  Summary 

Application  of  the  PAPA  detector  and  speckle  imaging  algorithms  at  CfA  have 
proceeded  from  a  development  phase  to  a  productive  scientific  program.  The  tools 
available  will  continue  to  improve,  both  on  the  analytic  side  and  in  the  capability  of  the 
hardware.  However,  our  current  program  has  produced  extensive  scientific  results  on  a 
range  of  astronomically  interesting  problems  including  supergiants  (Karovska,  et  al, 
1986a),  Young  Stellar  Objects  (Nisenson  et  al,  1985),  the  supernova  SN1987A 
(Nisenson  et  al,  1987,  Karovska  et  al,  1988),  Mu  Cas  (Karovska  et  al,  1986b)  and 
Seyfert  nuclei  (Ebstein  et  al,  1989).  The  crucial  component  appears  to  be  the 
availability  of  a  good  detector  that  is  photon  noise  limited,  for  which  the  photon  noise 
biases  in  the  reconstruction  processes  can  be  accurately  corrected. 


'll,  .  author  wishes  to  acknowledge  the  many  contributions  made  to  the  speckle 
imaging  program  at  CfA  by  Costas  Papaliolios,  Steven  Ebstein,  Margarita  Karovrka, 
Robert  Stachnik,  Robert  Noyes,  Clive  Standley,  Laurent  Koechlin,  Nathaniel  Carleton, 
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